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SURVEY 

Even  though  tins  is  a short  Progress  Report,  there  are  a number  of  interest- 
ing results  which  have  appeared  during  the  period  which  it  covers.  One  ol  these  is 
the  calculation  of  the  density  of  states  curve  for  nickel,  which  Hosier  has  been  carry- 
ing out.  Much  use  had  been  made  in  the  past  of  the  density  of  states  curve  for  this 
element  which  Krutter  and  1 calculated  many  years  ago.  This  curve  showed  a dip  in 
the  middle  of  the  3d  band,  which  w s of  considerable  importance  in  the  theory  of  ferro- 
magnetism, electronic  specific  heat,  and  other  properties  of  the  transition  elements. 
More  recent  work,  both  by  the  cellular  and  the  tight  binding  approximation,  had  shown 
that  the  details  of  the  old  cellular  calculation  on  which  this  work  of  Krutter's  and  mine 
was  based,  were  certainly  not  correct:  the  detailed  curves  of  energy  versus  k,  which 
now  appear  to  be  fairly  certain  from  newer  calculations,  arc  quite  different  from  our 
older  ones,  as  a result  of  the  small  number  of  spherical  harmonics  which  were  used 
in  this  older  cellular  calculation.  As  mentioned  in  the  preceding  Progress  Report,  it 
seemed  worthwhile  to  make  a new  calculation  of  density  of  states  for  nickel,  on  the 
basis  of  the  tight  binding  calculation  of  Fletcher  and  VVohlfarth,  which  seems  to  agree 
qualitatively  rather  well  with  the  new  calculations  on  copper  made  by  Howarth,  by  both 
the  improved  cellular  method  and  the  augmented  plane  wave  method.  Fletcher  and 
Wohlfarth  computed  the  density  of  states  just  at  the  top  of  the  energy  band,  but  the 
mathematical  problem  was  too  severe  for  them  to  wish  to  carry  out  ihe  complete  cal- 
culation. 

The  facilities  of  Whirlwind  make  this  calculation  rather  easy,  and  consequently 
Foster  has  carried  through  this  work,  and  reports  on  it  in  the  present  Progress  Re- 
port. His  result  is  very  interesting,  in  that  it  shows  a curve  closely  resembling  the 
old  one  which  Krutter  and  1 had  found.  The  dip  in  the  center,  which  was  the  basis  of 
most  of  the  qualitative  deductions  from  the  curve,  remains  much  as  in  the  earlier  cal- 
culations. Hence  we  may  consider  these  many  deductions  from  the  earlier  curve  to 
be  still  valid,  even  though  the  details  of  the  earlier  calculations  were  not.  What  has 
happened  is  that  the  errors  in  the  earlier  calculation  become  rather  well  ironed  out 
in  the  process  of  taking  the  density  of  states  curve:  we  have  about  as  many  states  as 
before  in  each  energy  range,  even  though  they  are  located  in  different  parts  of  the 
Briilouin  zone. 

I have  just  mentioned  the  calculations  of  Howarth.  He  left  to  return  to  Eng- 
land before  being  able  to  finish  completely  his  work  on  copper,  but  he  obtained  the 
essential  results  required  from  Whirlwind,  and  will  shortly  be  able  to  finish  his  cal- 
culations. The  net  result  was  that  he  found  that  the  augmented  plane  wave  method 
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converges  remarkably  rapidly,  more  so  than  the  cellular  method,  and  it  seems  as  of 
n<.yvv  to  be  the  best  meihod  which  we  have  for  solving  the  periodic  potential  problem. 

1 he  qualitative  results  on  copper  resemble  those  of  Fletcher  and  Wohlfarth  for  nickel, 
as  determined  by  the  tight  binding  method.  The  agreement  between  the  results  for 
copper  by  the  augmented  plane  wave  method,  and  by  the  cellular  method,  is  not  per- 
fect. This  may  be  partly  a result  of  the  fact  that  the  calculations  were  not  made  for 
exactly  the  same  potential,  for  evidence  is  rapidly  increasing  that  energy  bands  are 
remarkably  sensitive  to  small  changes  in  potential.  To  check  this  point,  Howarth  is 
preparing  to  repeat  the  cellular  calculation,  for  exactly  the  same  potential  used  for 
the  augmented  plane  wave  method,  so  as  to  get  a precise  comparison  of  the  two  methods; 
he  finished  the  necessary  calculations  for  this  on  Whirlwind  before  he  left.  We  shall 
hope  to  learn  more  details  of  the  resu'ts  later. 

Another  calculation  which  has  led  to  a definite  result  is  that  of  Kaplan  on  the 
ammonia  molecule.  Just  before  he  left  to  take  a new  post  at  the  University  of  Buffalo, 
he  obtained  preliminary  results,  also  by  means  of  Whirlwind,  for  the  LCAO  SCF  cal- 
culation for  this  molecule,  on  which  he  has  been  working  for  a long  time;  though  there 
is  still  considerably  more  to  be  done  before  the  work  is  completed.  These  preliminary 
results  indicated  a binding  energy  about  75  percent  of  the  observed  value.  This  repre- 
sents one  of  the  most  complicated  molecules  which  have  been  attacked  by  this  method; 
the  calculation  was  entirely  non-empirical,  and  the  values  of  the  three-  and  four-center 
integrals  were  calculated,  rather  than  being  estimated. 

McWeeny  has  also  left,  to  return  to  England,  and  before  leaving  he  finished 
up  the  calculations  on  the  water  molecule,  which  he  describes  in  this  Progress  Report. 
His  aim  was  to  find  whether  a simple  configuration  interaction  method,  based  on  the 
VB  method,  could  lead  to  results  significantly  better  than  the  LCAO  SCF  method.  The 
results  are  somewhat  discouraging:  he  was  not  able  to  do  nearly  as  well  as  the  LCAO 
SCF  method,  as  it  was  applied  to  this  molecule  by  Ellison  and  Shull.  Taken  together 
with  his  earlier  results  on  benzene,  and  with  Kaplan's  success  with  ammonia,  it  begins 
to  be  clearer  than  before  that  the  LCAO  SCF  method  is  probably  much  better  fitted  to 
a discussion  of  molecular  energies  and  wave  functions  than  other  easily  usable  methods. 
It  is  fortunate,  in  this  connection,  that  Meckler's  Whirlwind  program  for  the  LCAO 
SCF  method  is  available.  This  program  is  proving  to  be  remarkably  flexible;  and 
Meckler,  in  the  present  Progress  Report,  outlines  its  uses,  and  those  of  some  of  the 
other  programs  which  are  now  being  used  on  Whirlwind. 

Lowdin,  in  the  last  preceding  Progress  Report,  outlined  several  interesting 
aspects  of  the  self-consistent  field  method;  more  extended  papers  on  the  same  subject, 
which  he  has  prepared  for  publication,  have  been  accepted  by  The  Physical  Review. 

In  connection  with  the  "natural  orbitals"  which  he  introduced  in  those  papers,  Koster 
examines  a simple  special  case  in  the  present  Progress  Report,  which  throws  some 
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doubt  on  the  practical  usefulness  of  these  natural  orbitals  in  obtaining  rapid  conver- 
gence in  a configuration  interaction  problem. 

A number  of  other  problems  under  consideration  are  coming  along  well. 
Schultz  is  getting  well  into  the  problem  of  the  interaction  of  electrons  and  lattice  vibra- 
tions, a very  involved  problem  which  will  probably  receive  much  more  attention  in  the 
future.  White's  work  on  the  elastic  vibrations  and  elastic  constants  of  copper  is  near- 
ing completion,  and  represents  a very  interesting  piece  of  work.  Wood’s  self-consist- 
ent calculation  of  the  iron  atom,  using  two  different  potentials  for  the  electrons  of  plus 
and  minus  spin,  is  practically  finished,  and  shows  that  the  potentials  and  charge  dis- 
tributions for  the  two  spins  are  not  as  different  as  one  might  have  thought.  Saffrer.  is 
coming  along  with  extensions  of  the  augmented  plane  wave  method,  and  Howland  is 
starting  consideration  of  the  application  of  this  method  to  a polar  crystal,  KC1.  Other 
work  is  also  progressing,  as  the  items  in  the  deport  will  indicate. 

1 have  already  staled  that  Kaplan,  tlowarth,  and  McWeeny  have  recently  left 
the  group,  after  profitable  periods  with  it.  White  expects  to  leave  shortly.  Pratt  and 
Kleiner  have  taken  up  their  duties  with  Group  35  at  the  Lincoln  Laboratory,  and  while 
we  see  them  frequently,  they  no  longer  belong  directly  to  this  group.  Schweinler, 
after  a very  profitable  summer  at  the  Oak  Kidge  National  Laboratory,  has  decided  to 
remain  there  permanently.  We  shall  shortly  add  to  the  group  K.  K.  Ncsbet,  who  took 
his  bachelor's  degree  at  Harvard,  and  is  just  finishing  his  work  for  the  doctorate  at 
the  University  of  Cambridge,  where  lie  has  been  working  with  Dr.  S.  F.  Boys  on  prob- 
lems very  similar  to  those  which  have  concerned  this  group.  During  the  summer,  we 
have  had  two  visitors:  Dr.  A.  Delgano,  of  the  University  of  Belfast,  who  has  been  at 
the  Institute  m connection  with  ti  e Foreign  Students'  Summer  Program;  and  Dr.  P. 
Merryman,  a member  of  Professor  Mulliken’s  group  at  the  University  of  Chicago, 
who  has  been  here  studying  the  application  of  Whirlwind  to  the  calculation  of  molecular 
integrals.  We  expect  a visit  shortly  from  Dr.  K.  Kuedenberg,  also  of  the  University 
of  Chicago,  also  in  connection  with  this  same  problem. 

J.  C.  Slater 
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1.  CALCULATION  OF  A DENSITY  OK  STATES  CURVE  FOR  NICKEL 


The  calculation  of  the  density  of  states  curve  for  nickel  which  was  described 
in  the  last  Progress  Report ^ has  now  been  completed.  The  density  of  states  curve 
is  shown,  in  Fig.  1-1.  The  curve  has  been  normalized  to  have  a total  area  of  five  cor- 


ing. 1-1 

Density  of  states  curve  for  nickel 

responding  to  five  states  of  a given  spin.  The  method  of  calculation  of  this  density  of 
states  curve  from  the  energies  of  the  d bands  calculated  on  the  10  degree  mesh  in  the 
first  Bnllouin  zone  is  given  in  a recent  paper  by  J.  C'.  Slater  and  the  author  on  the 
tight  binding  method.  ^ ' 

The  striking  feature  of  this  density  of  states  curve  is  the  dip  in  the  middle. 
This  curve  can  be  compared  with  the  density  of  states  curve  for  the  d bands  in  the 
body-centered  structure  recently  calculated^  using  the  same  parameters  for  the 
nearest  neighbor  interactions  used  in  this  calculation.  This  curve  also  showed  a pro- 
nounced dip  in  the  center  of  the  distribution.  Thus  it  seems  on  the  basis  of  this  sim- 
plified tight  binding  calculation  that  both  the  body -centered  and  the  face -centered 
structur  e show  a minimum  for  the  density  of  states  curves  at  some  point  near  the 
middle  of  the  d bands.  This  dip  in  the  density  of  states  curve  for  a face -centered 
structure  was  predicted  many  years  ago  by  Knitter^  by  using  the  cellular  method 
and  it  is  interesting  to  find  that  the  calculation  carried  out  by  the  tight  binding  method 
shows  a similar  dip. 

The  purpose  of  this  calculation  was  to  extend  the  calculation  of  the  density  of 

(4) 

states  curve  given  by  Fletcher  and  Wohlfarth  over  the  entire  range  of  energies  for 
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the  d bands.  We  can  therefore  compare  only  a small  portion  of  the  density  o"  states 
curve  with  that  calculated  already.  In  the  figure  we  show  a smooth  curve  drawn 
through  the  step-like  curve.  This  eurve  extends  over  the  range  of  energies  that  Flet- 
cher and  Wohlfarth  covered  in  their  density  of  states  curve.  This  compares  quite 
well  with  the  curve  published  by  the  aforementioned  authors  which  tends  to  give  us 
confidence  in  the  method  of  calculating  the  density  of  states  we  have  used. 
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2.  ENERGY  BANDS  IN  COPPER 


The-  calculation  of  the  electronic  eigenvalues  of  metallic  copper  by  the  aug- 
mented plane  wave  method  has  been  completed  and  an  account  of  it  will  shortly  be 
submitted  for  publication.  The  rapid  convergence  reported  previously'  ' has  been 
found  in  every  state  considered.  In  all  cases,  it  has  been  found  that  combining  two 
solutions  inside  a sphere,  both  joined  to  the  same  plane  wave  outside  the  sphere,  is 
sufficient  to  give  the  "best"  wave  function  possible  using  only  one  plane  wave.  In 
other  words,  two  branches  of  the  E versus  Eq  curves^  are  sufficient  to  give  accur- 
ate eigenvalues  for  the  lowest  states.  Inclusion  of  six  such  terms  gives  no  change  in 
the  lowest  eigenvalues,  and  yet  does  not  produce  satisfactory  convergence  in  higher 
eigenvalues.  The  method  would  therefore  appear  to  be  limited  at  present  to  calcula- 
tions of  states  in  the  valence  and  conduction  bands,  and  to  obtain  accuracy  in  these, 
the  numerical  labor  is  not  excessive. 

In  the  case  of  copper,  the  inclusion  of  plane  waves  corresponding  to  wave 
vectors  in  t lie  two  innermost  Brillouin  zones  is  sufficient  to  produce  convergence. 

This  is  not  surprising  since  (1)  the  conduction  electrons  of  copper  deviate  only  slightly 
from  the  free  electron  behavior  and  (2)  the  wave  functions  of  the  u electrons  are  al- 
most entirely  confined  inside  the  atomic  sphere. 

Some  of  the  results  obtained  for  copper  are  given  in  Tables  2-1  and  2-2  below;  * 


labia  2-1 

Eigenvalues  of  conduction  band  of  copper  using  (a)  Hartree  and 
(b)  Hartree -Foek  atomic  potentials 
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*The  notation  used  is  fully  described  in  Ref.  3.  Eigenvalues  are  quoted  in  Rydberg 
units. 
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others  are  being  analyzed  at  the  time  of  writing.  As  an  example  of  the  power  of  the 
method,  a point  of  no  symmetry  in  k-space  lias  been  considered,  the  point  A,  (— , — , 
Using  two  solutions  inside  the  sphere  for  each  plane  wave,  and  using  7 plane  waves 
(with  wave  vectors  in  the  two  innermost  zones)  a secular  equation  of  order  14  was 
solved  in  this  case,  the  highest  order  equation  which  has  occurred  in  this  work.  The 
power  of  the  method  may  be  judged  from  the  results;  in  all  cases,  the  error  due  to 
truncation  of  the  infinite  series  of  augmented  waves  is  estimated  as  less  than  0.  5 in 
the  last  figure  quoted. 

Comparison  of  these  results  with  those  obtained  for  copper  by  the  cellular 
(3) 

method  shows  considerable  discrepancies,  particularly  in  the  position  of  the  states 
lying  in  the  d-band.  The  total  width  of  the  tl-band  remains  approximately  the  same 
(4.  I c.  v.  in  the  present  work),  but  the  top  of  the  band  lies  at  the  center  of  the  Brillouin 
zone,  the  base  lying  at  the  end  of  the  (110)  axis.  The  reason  for  the  discrepancy  in 
the  two  sets  of  results  appears  to  lie  in  the  different  potentials  used;  a cellular  calcu- 
lation of  the  eigenvalues  at  the  center  of  the  zone  in  the  d-band  has  been  carried  out 
using  the  exact  potential  used  in  the  present  work  differing  from  that  used  in  the  ori- 
ginal cellular  calculation  in  that  it  is  a constant  outside  the  atomic  sphere.  The  agree- 
ment between  this  calculation  and  the  results  of  the  augmented  plane  wave  method  is 

(3  4) 

quite  satisfactory.  Hence,  as  has  oeen  observed  in  other  work,  ’ the  excited 
energy  levels  in  crystals  appear  to  be  extremely  sensitive  to  the  crystal  potential 
used,  and  to  obtain  results  of  real  physical  significance,  an  attempt  at  self-consistency 
would  seem  necessary. 

In  concluding  this  series  of  reports,  I would  like  to  express  my  deepest  grati- 
tude to  Professor  J.  C.  Slater  for  making  possible  my  stay  with  the  Solid-State  and 
Molecular  Theory  Group,  and  for  his  continued  advice  and  encouragement  throughout 
the  year.  My  thanks  are  due  to  the  whole  group  for  many  stimulating  discussions  and 
for  much  helpful  advice.  The  entire  numerical  work  involved  in  this  project  has  been 
carried  out  on  the  Whirlwind  I digital  computer,  and  I am  deeply  indebted  to  all  the 
staff  there  for  their  help  and  cooperation.  Availability  of  the  computer  at  this  time 
was  made  possible  by  the  Office  of  Naval  Research,  to  whom  thanks  are  due. 
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i.  AN  AUGMENTED  PLANE  WAVE  METHOD  A3  APPLIED  TO  SODIUM 


Howarth's  program  for  the  augmented  plane  wave  method  having  been  applied 
to  copper  is  now  being  applied  to  sodium.  In  the  meantime  I have  been  attempting  to 
modify  Howarth's  program  so  that  calculations  on  the  energy  bands  in  sod'iim  can  be- 
gin according  to  a modification  of  the  APW  method  which  I outlined  in  my  last  Report.  ^ 

The  modification  requires  only  that  the  energy  expression  used  by  Howarth  be  replaced 

i ’ \ 

by  the  somewhat  more  complex  expression  which  is  given  in  my  last  Report.  v‘ ' I have 
also  undertaken  to  build  upon  the  part  of  Howarth's  program  which  (in  essence)  makes 
properly  symmetrized  combinations  of  APW  before  the  APW  are  allowed  to  interact  in 
the  secular  equation;  the  plan  is  to  simplify  the  coding  of  symmetry  points  in  the 
Hrillouin  zone  so  as  to  make  more  complete  the  mechanization  of  the  steps  of  the 
method  which  utilize  the  operations  of  the  cubic  point  group. 

The  purpose  of  the  present  calculations  is  to  examine  the  energy  bands  of 
sodium  and  also  to  compare  the  convergence  of  the  APW  scheme  with  the  modified 
scheme. 
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4.  ENERGY  BANDS  IN  POTASSIUM  CHLORIDE 


It  is  proposed  to  calculate  the  energy  band  structure  of  potassium  chloride 
using  the  augmented  plane  wave  (APW)  method^'  to  obtain  the  energies  for  certain 
high  symmetry  values  of  the  propagation  constant  k,  and  using  the  Slater -Koster  tight 
binding  scheme^  for  interpolation.  A rigid  lattice  will  be  assumed,  and  on  this  as- 
sumption excited  levels  will  be  included.  The  works  of  S.  R.  Tibbs^  and  of  W. 

(4) 

Shockley  provide  some  preliminary  information  about  alkali  halide  energy  levels 
and  functions,  and  the  recent  work  of  P.  O.  L.owdin^  on  cohesive  energy  provides 
tools  which  will  be  quite  useful  in  the  proposed  band  calculation.  With  the  information 
resulting  from  such  a calculation  it  is  hoped  that  new  light  may  be  shed  on  such  prob- 
lems as  the  importance  of  lattice  interactions  in  optical  transitions  and  the  effect  of 
impurities  (such  as  the  K-ccntcr)  on  the  energy  band  structure.  There  is  no  present 
work  which  attempts  a complete  quantum -mechanical  calculation  of  the  desired  energy 
levels  and  functions. 

The  present  attempt  at  a complete  calculation  is  made  possible  by  the  rapid 
convergence  of  the  new  APW  method  of  Slater  and  Saffren.  The  value  of  the  method  is 
illustrated  by  the  work  of  D.  J.  Howarth^  on  the  energy  band  structure  of  copper. 
Certain  difficulties  arise  in  adapting  the  method  for  an  ionic  crystal,  however.  First, 
the  potential  which  is  to  be  used  must  be  spherical  within  a sphere  about  each  lattice 
site  and  constant  outside  the  spheres;  in  our  case  the  real  potential  will  probably  be 
quite  non-sphcrical  and  non-constant  in  these  regions.  Secondly,  the  computer  pro- 
gram for  Whirlwind  developed  by  llowarth  for  the  copper  calculations  must  be  suffi- 
ciently generalized  to  take  account  of  the  two  types  of  centers  in  the  alkali  halide  unit 
cell.  As  Howarth  has  already  shown,  the  generalization  of  the  first  part  of  his  pro- 
gram will  not  be  difficult;  that  part  consists  of  obtaining  the  E versus  Eq  curves  of 
the  APW  method  and  the  resultant  single  APW's.  On  the  other  hand,  the  second  part, 
solving  the  secular  equation  between  interacting  APW's,  will  probably  be  quite  diffi- 
cult due  to  the  fact  that  the  linear  combinations  of  APW's  will  not  be  real  as  they  were 
for  Howarth.  This  problem  is  presently  being  studied. 

The  first  difficulty  inherent  in  the  proposed  calculation  concerns  the  initial 
crystal  potential  which  is  to  be  used.  The  model  we  are  employing  is  that  of  a rigid 
lattice  with  "free  ions",  Cl  and  K + , located  at  alternate  lattiee  sites.  The  electron 
eharge  density  of  the  crystal  is  taken  as  that  which  would  exist  if  all  the  electrons 
were  located  in  Hartree-Hartree  determined^  free  ion  functions  centered  on  the  vari- 
ous lattice  sites. 

The  ground  state  + of  the  crystal  is  written  as  a single  determinant  of  elec- 
trons in  free  ion  states;  this  is  possible  because  the  ions  K+  and  Cl  both  tmve  closed- 
shell  elec*ronic  structures,  with  eighteen  electrons  apiece.  The  charge  density  p(r) 
is  obtained  by  integrating  - over  the  coordinates  of  all  but  one  electron  and  sum- 
ming over  all  spins.  The  usual  periodic  boundary  conditions  are  used,  the  periodic 
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cell  being  a cube  containing  N C'l  ions  and  N ions.  ifi  is  normalized  so  that  the 
integral  of  p(r  ) over  this  cube  is  - 36  Me. 

Written  out  to  the  first  order  in  overlaps,  the  charge  density  is: 


PU  , l 


Pfrec(rl)  + p 


overlap  I 


(*•,)  • 


(4.  1) 


Here  Pfl.ec^ril  's  that  charge  density  which  is  obtained  by  simply  summing  the  charge 
densities  - eu*u  at  r of  the  electrons  in  their  free  ion  functions  u,  with  an  appropri- 
ate normalization  factor;  ami  Povorj  ap(r  l ) IS  Proportional  to  the  following  sum  of 
products: 


y1  s(qm|q  'm')  f X 
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Here  u *(  I ) is  the  free  ion  function  with  quantum  numbers  m,  centered  on  lattice  site 
in  1 

q,  and  containing  electron  1,  s is  the  overlap  integral  (including  spin  summation),  and 


the  prime  on  the  summation  symbol  means  (qm)  f (q'm1).  All  the  terms  in  p 


involve  just  two  sites  and  contain  all  functions  u as  sums  of  the  form  Y u cl(  1 ) 


overlap 


(5)  m 

v 1 (which  follows  directly  from  the  additional  theorem 


u *(2).  Bv  a theorem  of  L.bwdin 
in 

for  spherical  harmonics)  each  such  term  may  be  evaluated  by  taking  the  vector  be- 
tween the  two  sites  as  a z-axis  for  the  two  ions  affected.  On  this  basis  the  overlap 
integrals  calculated  by  Lowdin  for  KC1  may  be  used  directly,  and  the  evaluation  of  p 
should  not  be  too  difficult. 

The  Hartree-Foek  equation  for  the  one-electron  function  q>(r)  can  be  written 
in  the  following  form: 
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h_ 
2 m 


e V . (r  ) - e V , (r  ) 
eoul  exch 
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(4.  2) 


The  potential  which  must  he  evaluated,  then,  is  V(r  ) = V(  j(r)  + V ^ 
exchange  potential  V ^ may  lie  approximated  by  use  of  Slater's  p 1 / ^ equation, 
where  p is  calculated  as  described  in  the  preceding  paragraph. 


(r).  The 
(8) 


V 


coul 


(r)  is  made  up  of  contributions  from  both  the  electron  cloud  and  the 


nuclei.  By  various  elementary  operations  ^couj(r)  can  written  as  follows  (in  the 
MKS  system): 


V ,(r  ) = 

COUl 


4tte  a 
o 


^S(r  ) + u'l(r  ) + uiF(r)j, 


(4.  3) 


where  a is  the  interionic  distance,  and  r is  measured  in  units  of  a.  The  terms  on  the 
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l ight  of  (4.  i)  are  defined  as  follows: 


S(r>  = £ 


±lf 

I r - H 


(4.  4) 


where  K is  a vector  to  the  lattice  site  q,  q even  being  a Ci  site,  and  q odd  being  a 

+ q 

K site; 
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(4.  5) 


where  U (x)  is  the  total  radial  charge  density  at  a radius  x in  a closed  shell  ion  of  the 
type  which  is  located  at  lattice  site  q; 
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dv,’ 


where  q = u is  ihc-  origin  (a  <JT  site)  and  q = 1 is  the  nearest  K+  site  in  the  -x  direc- 
tion; and 
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(4.  7) 


Physically,  S(r  ) (Eq.  4.  4)  is  the  potential  at  a point  r in  a cubic  lattice  of 
unit  point  charges  with  alternating  sign  if  one  of  the  negative  unit  charges  is  taken  as 
an  origin.  From  symmetry  it  can  be  shown  that  S(x,  y,  z)  = - S(  1 - x,  y,  z).  S(r  ) is 
presently  being  evaluated  in  the  central  unit  cell  by  expanding  S(i  ) - p in  a Taylor 

series  about  r = 0;  the  coefficients  of  the  products  xa  y^  zC  take  the  form  of  combina- 

(9) 

tions  of  lattice  sums.  The  latter  are  being  calculated  by  Evjen's  method.  This 
method  is  not  sufficiently  accurate  for  all  tlm  sums  involved,  however;  hence  a few 
of  the  sums  will  be  calculated  by  an  extension  of  the  Ewald  method.  S(r  ) may  also 
be  checked  by  direct  use  of  the  Ewald  method^ at  specific  points  in  the  lattice. 
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[)  q 

For  (a  ' b + c)  odd,  the  coefficients  of  >:  y z are  of  course  7ero;  it  appears  that 
satisfactory  accuracy  is  obtained  when  the  expansion  is  extended  to  include  the  tenth 
order  terms.  The  worst  point  in  the  unit  cell  is  (x,  y,  7.)  = (1  /2,  l/2,  1 /2),  but  this 
point  is  not  included  in  either  of  the  touching  ionic  spheres;  the  larger  sphere  (Cl  ) 
intersects  the  ( 1 , 1 , 1 ) direction  at  ( 1/3,  1 / 3 , 1 / 3 ) . 

T(r)  (Eq.  4.  5)  represents  the  correction  to  the  potential  which  arises  from 
the  spatial  extension  of  the  electronic  charge.  The  contributions  of  nearest  neighbors 
are  large,  but  those  of  the  further  neighbors  fall  off  rapidly  with  distance  from  the 
origin. 

F(r)(Eq.  4.6)  mostly  represents  the  potential  contribution  arising  from  the 
overlap  charge  density  (Eq.  4.  !).  Several  approximations  for  this  term  are  presently 
being  studied,  but  the  term  has  not  yet  been  evaluated.  The  terms  in  F(r  ) are  all 
such  that  the  Lowdin  theorem  is  applicable,  and  the  Lowdin  overlap  integrals  can 
again  be  used  directly. 

u(Eq.  4.  7)  is  the  charge  normalization  factor;  it  is  easily  evaluated  from 
Lowdin's  overlap  integrals.  Of  the  four  terms  contributing  to  the  potential,  only  S(r  ) 
appears  without  a factor  w. 

Treating  the  component  terms  as  indicated  above,  V(r)  will  be  evaluated  on 
the  surfaces  of  spheres  centered  on  lattice  sites  with  radii  equal  to  the  empirical 
ionic  radii.  The  deviation  from  sphericity  will  be  studied.  If  it  does  not  seem  too  un- 
reasonable, the  potential  V(r  ) will  be  replaced  by  its  spherical  averages  within  the 
spheres  and  by  an  effective  constant  value  outside  them.  The  APW  method  will  be 
applied  using  this  as  an  initial  potential,  and  approximations  to  E and  <$>  (Eq.  4.  2) 
will  be  sought. 

KC1  has  been  chosen  as  the  particular  alkali  halide  to  be  studied  because  for 
it  the  free  ion  functions^  and  the  overlap  integrals^  have  already  been  calculated, 
and  because  the  ions  Cl  and  K+  are  of  comparable  size,  thus  putting  less  strain  on 
the  accuracy  of  the  expansion  of  S(r  ).  More  overlap  integrals  are  required  than  for 
NaCl,  for  instance,  but  the  overlaps  are  also  smaller.  These  reasons  do  not  prohibit 
calculation  for  the  other  alkali  halides,  however,  and  future  developments  may  even 
show  that  calculations  for  one  of  them  might  be  preferable. 

Finally,  in  anticipation  of  the  results  from  the  APW  method,  the  tight  binding 
interpolation  formulas  have  been  developed.  3s,  3p,  and  4s  orthogonal,  localized 
functions  were  assumed  for  each  type  of  ion.  The  4s  functions  were  included  to  allow 
consideration  of  excited  levels.  The  matrix  components  of  the  Hamiltonian  between 
Bloch  functions  were  written  out  with  the  matrix  components  between  localized  func- 
tions appearing  as  disposable  constants.  With  second-nearest-neighbor  interactions 
considered,  there  are  thirty-four  such  constants,  while  with  only  nearest  neighbors, 
there  are  eighteen.  Some  of  these  may  be  neglected,  however,  on  the  assumption 


A 
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that  the  orthogonalized  functions  are  about  the  same  as  the  free  ion  functions  they  re- 
place. I he  secular  equation  has  been  factored  at  many  symmetry  points  in  momentum 
space,  and  the  roots  have  been  extracted.  When  the  APW  energies  are  available  at 
these  points,  the  disposable  constants  will  be  determined.  The  exact  numbers  of  points 
and  constants  actually  used  will  of  course  depend  on  how  much  difficulty  is  encountered 
in  carrying  out  the  APW  method  for  KC1. 
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5.  ENERGY  RANDS  IN  GRAPHITE 


A tight  binding  calculation  has  been  made  on  Graphite  (described  in  the  pre- 
vious Report),  but  the  results  were  unsatisfactory.  Specifically,  for  certain  values 
of  the  reciprocal  lattice  vector  k,  the  overlap  matrix,  which  should  be  positive  definite, 
was  found  to  have  negative  eigenvalues.  It  is  believed  that  this  contradiction  arose  be- 
cause of  the  following  situation.  In  constructing  the  matrix  elements  S^,  the  terms 
arising  from  overlap  integrals  beyond  the  third  nearest  neighbor  distance  were  neg- 
lected; however,  the  overlap  integrals  used  were  found  by  integrating  analytically 
over  all  space.  Thus,  failure  to  include  enough  neighbors  in  the  evaluation  of  the 
overlap  matrix  elements  allows  the  possibility  of  negative  overlap  eigenvalues,  an 
effect  analogous  to  an  insufficient  number  of  terms  in  the  Fourier  series  of  a positive 
definite  quantity. 

At  present,  the  calculations  are  being  re-examined  with  the  view  of  increas- 
ing the  number  of  neighbors  considered. 


F.  J.  Corbato 
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6.  A REVIEW  OF  T11E  USE  OF  WHIRLWIND  BY  THE  SOLID-STATE  AND 

MOLECULAR  THEORY  GROUP 

This  report  will  summarize  the  use  of  Whirlwind,  the  computing  machine, 
by  members  of  this  group.  Although  special  programs  have  been  written  for  the  mach- 
ine, programs  which  perform  a particular,  formalized  set  of  operations  specific  to 
one  problem,  the  emphasis  here  will  lie  on  the  description  of  more  general  programs: 
either  basic  mathematical  routines  which  are  called  for  naturally  in  the  calculus  of  a 
class  of  physical  theories,  or  programmed  physical  theories  which  are  standards 
ready  to  compute  on  any  structure,  needing  only  the  specification  of  an  atomic  number, 
potential,  or  a table  of  integrals.  Whirlwind  does  have  features  that  other  machines 
do  not,  but  no  program  has  been  written  for  Whirlwind  only;  except  that  with  the  mach- 
ine's high  speed  and  great  capacity,  we  may  be  singularly  adjusted  to  such  things  as 
large  order  matrices,  Highly  iterative  schemes,  and  wide  gamuts  of  parameters. 

The  Computer 

Whirlwind  1 is  a high-speed  electronic  digital  computed  located  at  M.  I.  T.  and 
sponsored  by  the  O.  N.  R.  It  was  built  and  is  maintained  as  a prototype  for  engineers 
and  mathematicians  in  the  development  of  circuitry  and  logic  for  digital  computers. 
Time  has  been  made  available  by  the  M.  l.T.  Committee  on  Machine  Methods  of  Com- 
putation and  Numerical  Analysis  (Professor  P.  M.  Morse,  Cha.rman),  which  super- 
vises the  scientific  and  engineering  applications  of  the  machine,  to  staff  membei  s and 
students  at  M.  1.  T. 

The  direct,  or  fast  storage  elements  of  the  machine  are  magnetic  cores. 

These  are  grouped  into  registers  each  of  which  can  accommodate  a word  length  of  16 
binary  digits.  There  are  2048  magnetic  core  registers  and  it  is  in  these  that  sections 
of  the  program  in  operation  are  stored  along  with  numbers  to  be  operated  upon.  That 
is,  control  and  arithmetic  is  done  only  with  the  core  memory.  As  passive  storage, 
whose  contents  when  needed  must  be  transferred  into  the  cores  at  a relatively  slow 
rate,  there  are  reels  of  magnetic  tape  and  a 22,  528  register  magnetic  drum.  Another 
drum  will  be  available  soon  and  there  is  the  outside  storage  of  punched  paper  tape. 
Information  is  fed  into  the  machine  on  punched  paper  tape  passed  through  a photoelec- 
tric reader  and  the  output  can  be  had  on  automatic  typewriter,  punched  tape,  magnetic 
tape,  or  a photographed  cathode-ray  tube  display. 

The  electronic  speed  of  the  machine  is  about  40,  000  operations  per  second, 
per  se.  The  bare  electronic  machine,  however,  is  not  the  easiest  to  code  and  needs 
to  be  polysyllabic  in  its  words.  The  basic  word  length  corresponds  to  about  4 decimal 
digits  and  necessitates  fixed  decimal  (or  binary)  point  arithmetic.  Therefore,  a pro- 
grammed arithmetic  has  been  written:  numbers  cccupy  two  registers  in  the  form  x2^ 
and  the  order  to  multiply,  for  example,  transfers  control  to  a slew  of  operations  which 
perform  the  double-length,  floating  point  multiplication.  The  programmed  arithmetic 
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occupies  a section  of  fast  storage,  numbers  take  up  twice  the  space,  and  more  time  is 
spent  on  an  arithmetic  operation.  In  fact,  programmed  arithmetic  is  about  40  times 
slower  than  ordinary  Whirlwind. 

Associated  with  the  programmed  arithmetic  is  a wonderful  system  of  service 
routines.  These  include: 

Automatic  conversion  of  a program  written  in  almost  basic 
English  to  the  binary  code  of  the  machine. 

The  floating  address  system  whereby  a program  is  written 
with  addresses  referred  to  only  as  letters,  algebraic  symbols  whose 
values  are  assigned  by  the  machine  itself. 

Automatic  cycle  counters. 

A set  of  buffer  registers,  each  three  ordinary  registers 
long,  to  be  used  in  the  accumulation  of  partial  sums  without  the 
intermediate  round-off  into  two  registers. 

An  elaborate  system  of  program  checks  and  post-mortem 
routines  for  mistake  diagnosis. 

All'this  and  more  leaves  the  machine  with  sufficient  speed  and  capacity,  but  with  a 
logic  more  malleable. 

The  service  routines  are  developed  by  the  staff  of  mathematicians  at  Whirl- 
wind. Programs  for  calculation  are  written  by  the  originators  of  the  problems  who 
learn  coding  from  the  computer  staff  and  who  turn  to  the  staff  when  the  malicious  mach- 
ine reacts  absurdly  to  the  unique  and  careful  logic  of  their  programs.  The  staff  is 
most  capable.  It  is  always  shown  that  what  seems  a fundamental  incompatability  is 
a simple  lack  of  understanding  which  can  be  corrected. 

There  are  two  types  of  programs:  the  sub-routine  and  the  production  style. 

The  sub-routine  is  a packaged  unit  ready  to  be  inserted,  to  perform  a complex  opera- 
tion, into  any  large  program.  The  external  program  sets  up  the  data  as  required  by 
the  sub-routine,  transfers  control  to  the  sub-routine  which  does  what  it  should  and  re- 
turns control  to  the  external  program,  which  goes  on  from  there.  The  programmed 
arithmetic  is  a sub-routine,  and  many  of  the  matrical  routines  to  be  described  are  in 
sub-  outine  form,  easily  assembled  into  compound  programs.  The  other  type,  a pro- 
duction style  routine,  is  complete.  As  the  data  is  read  in,  the  program  assembles 
it,  has  it  manipulated  through  a series  of  smaller  sub-routines,  and  then  displays  the 
output.  Production  style  routines  are  not  meant  to  be  synthesized  into  larger  routines; 
sub-routines  are. 

The  AlgeDraic  Routines 

Much  of  matrix  algebra  has  been  programmed.  To  begin  with.,  there  is  the 
diagonalization  of  Hermitian  matrices.  This  method  is  this:  Let  H be  a real  sym- 
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* 


metric  matrix.  It  is  transformed  to  diagonal  form  by  a succession  of  2 x 2 orthogonal 
transformations.  Let  u be  an  orthogonal  matrix  with  unity  for  every  diagonal  element 
except  the  i.1*1  and  j1*1,  and  zero  for  every  off-diagonal  element  except  the  i.  and  the 
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c and  s are  determined  oy  the  trigonometric  identity 


2 

c + 


2 

s 


1 


and  the  condition  that  IF  = 0. 
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squares  of  all  the  elements  of 


Under  any  orthogonal  transformation  the  sum  of  the 
II  is  unchanged: 
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which  is  invariant.  However,  this  particular  transformation  has  increased  the  sum 
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of  the  squares  of  the  diagonal  elements  by  211..: 
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The  invariance  of  the  truer  and  determinant  imply 


Tr  H'  = Tr  H 


If.  If.  = H H . - H2 
•i  JJ  n JJ  >J 


sci  that 


Y (If  )2  = (Tr  H)2 
^ mm 
m 


I 11  11+  2H2 

mm  nn  ij 


m ^ n 
2 


= y (H  r + 2H 

^ ' mm'  i.j 


Weight  lias  been  shiftel  to  the  diagonal  and,  if  for  each  2X2  transformation  the  lar- 
gest off-diagonal  element  is  made  to  vanish,  the  process  will  converge.  The  total  or 
thogonal  transformation  is  the  product  of  all  the  2x2  transformations,  and  is  ealeu 
luted  in  parallel  with  the  diagonalization.  After  each  reduction,  and  beginning  with  a 
unit  matrix,  u is  multiplied  into  the  orthogonal  matrix  constructed  up  to  that  point, 
affecting  only  the  i'*'  and  j1*1  columns.  H is  considered  diagonal  when  the  absolute  lar 
gest  off -diagonal  element  is  less  than  a prescribed  criterion. 

In  ‘lie  machine,  the  upper  half  of  the  symmetric  matrix  is  stored  in  fast  stcr 
age,  the  orthogonal  matrix  is  kept  on  the  drum.  The  routine  is  amazingly  fast  and  ac- 
curate, and  exists  as  both  a sub-routine  and  a production  style  program. 

Complex  llermitian  matrices  have  to  be  handled  by  enlargement.  Let 

H = A + iB 


where  A is  real  symmetric,  B real  antisymmetric.  The  eigenvalue- eigenvector 
equation  reads  as 


(A  + iB)(x  + ly)  = \(x  + iy) 
or 


Ax  - By  = \x 
Bx  + Ay  = \y 

The  last  two  equations  can  be  written  in  enlarged  matrix  form  as 


A -B 

X 

X 

= K 

B A 

y 

y 
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so  that  a real  symmetric  matrix  is  to  be  diagonalized. 

The  routine  has  been  called  for,  of  course,  in  those  theories  which  approxi- 
mate a wave  function  as  a linear  combination  of  orthogonal  functions.  The  configura- 
tion interaction  in  the  oxygen  molecule,  considered  by  the  writer,  was  analyzed  on 
Whirlwind.  Kaplan's  treatment  of  the  ammonia  molecule  will  require  the  solution  of 
a 14  x 14  secular  equation.  For  his  tight  binding  computation,  Kostcr  took  the  sub- 
routine, added  a program  which  calculated  the  matrix  elements,  and  ran  through  220 
independent  points  in  k- space  for  the  body  centered  structure,  680  points  for  the  face- 
centered.  It  should  be  realized  that  the  ability  to  handle  large  matrices,  once  they 
exist,  is  not  enough.  There  is  the  awful  problem  of  the  matrix  elements,  as  to  com- 
putation and  transcription.  The  program  written  for  Schultz,  to  be  described  in  some 
detail,  is  lovely  in  the  absolute  minimum  of  input  required  for  the  mechanization  of 
the  theory. * 

The  existence  of  a fast  and  accurate  diagonalization  sub-routine  makes  feas- 
ible many  other  matrix  computations.  Determinants  can  be  evaluated  as  the  product 
of  the  eigenvalues.  It  seems  best  to  evaluate  any  invariant  of  a matrix  by  first  diag- 
onalizing the  matrix  and  then  forming  appropriate  products  of  the  eigenvalues.  The 
calculation  of  the  inverse  of  a matrix  (or  the  square  root  of  the  inverse)  has  been  pro- 
grammed in  this  way.  The  matrix  is  first  brought  to  diagonal  form,  the  reciprocals 
of  the  diagonal  elements  are  taken  (or  the  square  roots  of  the  reciprocals),  and  then 
this  diagonal  inverse  is  brought  back  to  non-diagonal  form  by  the  reverse  of  the  or- 
thogonal transformation.  This  inversion  routine  is  in  sub-routine  form  and  in  com- 
plete production  style  form.  As  a sub-routine  it  was  incorporated  into  the  impurity 
level  calculation  of  Koster.  He  recognized  a Green's  function  as  the  inverse  matrix 
(H  - E)  , where  H is  the  unperturbed  Hamiltonian  and  E is  the  sought  energy  level. 
The  problem  is  solved  in  the  space  of  a finite  number  of  Wannier  functions  so  that  H 
becomes  a finite  matrix.  For  the  body-centered  structure  16  values  oi  E were  in- 
serted for  each  of  the  220  points  in  k-sapee. 

The  square  root  inversion  routine  was  useful  in  the  orthogonalization  of  a set 
of  functions,  and  in  the  solution  of  the  secular  equation  arising  among  a set  of  non- 
orthogonal  functions.  If  a trial  wave  function  be  expressed  as  a linear  combination  of 
non-orthogonal  functions 


u.(x)  = £ v (x)  a^ 

j J 


*The  work  of  the  Whirlwind  users  mentioned  here  is  reported  in  the  Progress  Reports 
of  this  group  and  in  the  publications  listed  at  the  beginning  of  these  reports 
biarinually. 
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the  eigenvalue-eigenvector  equation  reads  as 


I H..  a..  = I S.  a,  E. 
j ‘J  jk  j lJ  Jk  k 


with  the  side  condition 


Y a*  S a = 6. 
mi  — 


m,  n 


mn  nj 


ij 


where  H..  is  a matrix  element  of  the  Hamiltonian  and  S . is  an  overlap  matrix  ele- 
>J  tj  K 

ment 


H.  = f v.*(x)  H v.(x)  dx 
ij  J 1 ' ' j' 

S. . = f v.*(x)  v (x)  dx 

ij  J 1 j 

In  matrix  notation  the  simultaneous  equations  are 

HA  = SAE 
AfSA  = 1 


Following  a suggestion  of  L.owdin  (M.  I.  T.  lectures)  define 


U 


S+  l!Z  A 


so  that 


U*U  = 1 

(S'  l!Z  H S'  1/,a)  U = UE 


and  an  ordinary  secular  equation  is  to  be  solved.  * A program  which  carries  out  this 
Lowdin  scheme  has  been  written.  As  a production  style  routine  it  was  used  by  McWeeny 
in  his  valence  bond  studies,  and  as  a sub-routine  it  was  used  by  Howarth  in  the  last 
stage  of  the  augmented  plane  wave  method,  and  by  Corbato  in  his  tight  binding  calcula- 
tion on  graphite. 


*S 


' / 2 t 

*'  is  defineable  since  S = V V where  V . = v.(x).  S is  therefore  Hermitian  and, 

xj  J 

if  the  columns  of  V are  independent,  positive  definite. 
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Integrals  that  are  needed  in  the  molecular  calculations  done  here  are  first 
evaluated  among  the  basic  set  of  atomic  orbitals.  The  molecular  orbitals,  out  of 
which  are  compounded  the  rnany-electron  functions,  are  linear  combinations  of  the 
atomic  orbitals.  The  integrals  which  appear  in  matrix  elements  are  integrals  among 
the  molecular  orbitals.  These  are  derived  from  the  atomic  orbital  integrals  by  means 
of  a congruent  transformation  U*  A U,  where  A is  a suitably  ordered  matrix  of  atomic 
integrals,  and  U is  a matrix  whose  elements  are  quadratic  forms  in  the  coefficients 
of  the  linear  combinations.  ^ A congruent  transformation  program  exists  in  produc- 
tion style;  it  has  been  used  by  McWeenj  and  by  Kaplan.  As  a prerequisite  to  it,  a 
matrix  multiplication  sub-routine  was  written. 

Special  Programs 

A colossal  program  to  carry  out  the  augmented  plane  wave  method  for  copper 
was  written  by  Howarth,  and  is  being  further  developed  by  Saffren  for  application  to 
sodium.  From  this  specific  program  certain  sub-routines  can  be  pulled  out  for  use 
elsewhere;  Simpson's  rule,  generation  of  P^(cos  9),  generation  of  spherical  Bessel 
functions,  and  Gauss -Jackson  numerical  integration. 

The  adaptation  to  Whirlwind  of  Roothaan's  self-consistent  L.C.  A.O.  method 
is  another  large  and  specific  program.  It  has  been  used  by  Kaplan  for  the  ammonia 
molecule  and  will  be  used  by  Allen  in  his  polarization  study. 

Finally,  there  is  the  program  for  the  electron-lattice  interaction  model  by 
Schultz.  The  model  treats  one  otherwise  free  electron  in  interaction  with  the  field  of 
two  degenerate  lattice  oscillators.  The  wave  function  of  the  system  is  expanded  in  a 
set  of  functions  each  belonging  to  the  same  total  momentum,  which  is  a good  quantum 
number,  and  each  distinguished  by  the  number  of  quanta  in  the  oscillator  field.  There 
are  three  parameters  that  Schultz  wishes  to  vary  independently  and  widely:  the  total 
momentum  (t|),  the  oscillator  momentum  (a),  the  coupling  constant  (y).  For  each  set 
of  q , a , y the  convergence  of  the  energy  levels  is  to  be  followed  as  states  of  higher 
occupation  numbers  are  brought  in.  The  energy  matrices,  though  large  and  many, 
can  be  computed  easily  hy  hand;  there  are  an  enormous  number  of  zero  elements  and 
not  many  independent  varieties  of  the  non-zero  ones.  The  matrix  which  allows  up  to 
three  quanta  in  the  field  is  shown  in  Fig.  6-1.  The  angle  shaped  lines  are  boundaries 
between  states  of  different  occupation  numbers.  Each  occupation  number  brings  in 
that  number  plus  one  new  states  --  the  number  of  ways  to  distribute  the  quanta  be- 
tween the  two  oscillators. 

Now,  all  these  zeros,  which  present  no  computational  problem,  nevertheless 
would  have  to  be  written  down  on  tape  preparation  forms,  typed  and  punched  out  on 
paper  tape,  passed  through  the  photoelectric  reader  on  equal  footing  with  the  other 
numbers.  Easy  arithmetic  though  it  may  involve,  the  hand  computation,  transcription, 
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typing,  and  read  in  of  the  matrices,  would  be  time  consuming  and  greatly  susceptible 
to  error.  The  matrices  must  be  formed  in  the  machine.  Rut  is  the  machine  to  go  along 
through  the  matrix,  set  itself  at  an  i,  j position,  ask  if  the  number  to  be  inserted  is 
zero,  and  if  it  is  not,  is  the  machine  then  to  refer  to  a dictionary,  somehow  planted  in 
storage  which  will  direct  it  to  a set  of  operations  suuabie  lo  that  value  of  i,  j?  This 
now  rhetorical  question  was  answered  with  a no.  A block  of  registers  can  be  cleared 
immediately  by  the  program.  There  is  no  need  to  insert  zero  elements  or  even  to  pass 
through  ihem.  Rather  than  pick  a position  and  deeide  what  is  to  be  plaeed  there,  pick 
up  what  is  to  be  placed  and  decide  where.  This  is  how  the  program  works:  The  prob- 
lem has  been  solved  up  to  a certain  number  of  quanta  and  now  the  matrix  is  to  be  aug- 
mented by  the  introduction  of  states  with  one  more  quantum.  Let  m be  this  number 
of  quanta  plus  one.  The  order  of  the  matrix 

m(m  + 1 ) 

n = — i L 

Z 

The  new  off-diagonal  elements  are  all  of  the  form  \/ k y 


where 


k 


max 


m 


1 


and 


0 sf  q m + 2 


Each  k goes  into  two  positions  given  by 


^q  ^o  ^ 

V = j0  - (m  - 1)  + q 


i =1  - m 

q Jq 

i ' = j ' - (m  + 1) 

q Jq 


where 


i = n - 1 
Jo 


The  gross  eyele  is  the  one  on  m as  it  is  set  to  increasing  values  beginning 
with  m = 2.  (The  case  of  no  quanta  is  a special  one  whose  1 x 1 matrix  is  computed 
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Fig.  6-1 

Schultz's  matrix  for 
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immediately.  ) Rick  an  in.  From  this  find  n and  k . Establish  i . Set  up  a q 

max  Jo  r m 

counter  and  place  the  off-diagonal  elements.  The  diagonal  elements  are  of  the  form 

t * [l  + (^  ■ 2q)aJ2 

where  t = 2(m  - 1)  and  0 $ q ^ ^ . The  address  of  the  diagonal  element  is  n - 1 - q, 
n - 1 - q. 

The  program  input  consists  of  m , q,  a,  y.  m begins  at  2 and  runs  up  to 

the  m decided  on  by  Schultz.  Because  the  triangular  Hermitian  matrix  has  to  be 
max  J b 

in  fast  storage,  and  because  it  would  be  costly  timewise  to  break  up  the  program  into 

small  bits  which  are  transferred  from  the  drum,  the  m cannot  exceed  7,  which 

max 

corresponds  to  a maximum  matrix  order  of  28. 

The  only  apology  that  1 can  offer  for  going  into  the  pattern  of  this  program  is 
that  il  was  exciting  to  analyze,  and  I hope  that  the  inversion  of  thought  in  the  matrix 
formation  will  prove  to  be  good  experience  for  other  problems. 

Reference 

1.  This  technique  is  described  in  the  Quarterly  Progress  Report,  Solid-State  and 
Molecular  Theory  Group,  M.I.T.,  April  15,  1954,  p.  28  and  p.  32. 

A.  Meckler 
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SOLID-STATE  PHYSICS 

In  the  following  report,  various  applications  of  the  Whirlwind  1 computer  to 
problems  of  molecular  and  solid-state  physics  are  discussed  in  an  attempt  to  illustrate 
the  high  versatility  of  such  a machine.  Inasmuch  as  the  potentialities  of  high-speed 
computers  are  still  not  too  well  defined,  the  numerical  techniques  are  sketched  in 
some  detail  in  order  to  indicate  the  scope  of  reasonable  coding  procedures.  It  is 
hoped  that  in  this  indirect  fashion,  it  will  be  possible  to  display  some  of  the  shifts  in 
numerical  emphasis  allowed  by  a modern  computer. 

Function  Generators 

In  coding  a problem  for  a digital  computer,  there  often  arises  the  need  of 
knowing  the  value  of  a function  for  a given  argument.  There  are  two  usual  solutions 
to  this  problem:  One  may  either  use  a table  of  the  function  values  and  an  interpolation 
scheme  if  the  values  of  the  argument  are  non -r.nifor in,  or  one  may  use  a generation 
subroutine  each  time  to  form  the  function  value  from  the  argument.  If  the  function 
varies  rapidly  over  the  range  of  the  argument,  or  if  many  functions  are  required,  the 
table  method  will  require  a great  deal  of  the  storage  space  in  the  computer.  On  the 
other  hand,  generation  of  a complicated  function  can  be  a time-consuming  process  even 
on  a high-speed  computer.  In  any  case  the  problem  of  reading  tables  of  function 
values  into  the  computer  is  a serious  one,  for  the  unsatisfying  task  of  proofreading 
is  the  only  easy  way  to  eliminate  errors  among  the  entries.  Thus  function  generators 
are  extremely  useful  in  the  coding  of  most  problems,  either  in  direct  use  or  as  the 
preliminary  generators  of  tables  for  interpolation. 

Because  of  their  widespread  application,  the  ordinary  functions  such  as  ^/x", 
sin  x,  cos  x,  sinh  x,  cosh  x,  exp  x,  and  tn  x are  available  in  the  Whirlwind  Subroutine 
Library  which  has  been  prepared  by  the  computer  staff.  In  general,  when  coding 
problems,  these  relatively  simple  functions  are  best  obtained  by  the  direct  use  of  the 
generation  subroutines. 

In  the  writer's  current  work  on  the  energy  bands  of  graphite,  the  need  arose 
for  evaluating  various  types  of  two-center  integrals  involving  analytic  Slater  AO's. 

By  the  use  of  prolate  spheroidal  coordinates,  these  integrals  were  all  simple  ex- 
pressed in  terms  of  the  two  sets  of  functions 
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n:  0(1  )N 


0 < x < oo 
- oo  < y < oo 


vhere  x,  y,  and  N are  fixed  for  any  particular  integral,  Since  all  ordinary  two-center  in- 
tegrals of  the  Slater  AO’s  (except  exchange  and  hybrid  integrals)  are  simply  expressed  in 
terms  of  the  Cn  and  Bfil  a general  subroutine  has  been  programmed  for  the  generation 
of  the  two  sets  of  functions  for  given  x,  y and  N.  Extensive  application  of  this  sub- 
routine has  been  made  by  the  writer  and  by  Dr.  A.  Dalgarno  of  Queen’s  University 
Northern  Ireland  who  was  visiting  M.  l.T.  this  past  summer. 

It  is  perhaps  interesting  to  examine  in  some  detail  the  method  used  in  coding 
the  above  subroutine,  since  it  is  illustrative  of  the  numerical  techniques  which  are 
applicable  in  a high-speed  computer.  The  functions  Cn(x),  being  simply  polynomials 
in  x which  is  always  positive,  are  sums  with  all  positive  terms,  and  thus  offer  no  nu- 
merical cancellation  problems.  For  greater  computational  efficiency  though  it  is 
convenient  to  use  the  recursion  formula 


Cn(x)  = xn  + n Cn  _ t(x) 


to  form  the  functions  of  higher  n successively  from  CQ(x)  = 1.  The  functions,  Bn(y), 
however,  are  considerably  more  difficult  in  that  for  the  magnitude  of  y approximately 
less  than  n,  the  closed  form  given  above  has  large  numerical  cancellations.  To 
avoid  this  the  following  reformulation  was  made: 


Bn(y)  = 2(-  1)“ 


n 

«2(-l)r,ni  I* 
1 = 0, 


(2j  * 1) 

(n  - j)!  ! (n  + j + 1).’ 


‘iw 


where 


i (x)  s (-  t,)n  j (tx)  = spherical  Bessel  function  of  the  first  kind  with 

imaginary  argument  (i.  "modified  spherical 
Bessel  function") 

n!  .'  = n(n  - 2)(n  - 4)..,  (2  or  1) 

Y'  = sum  over:  j even  if  n even;  j odd  if  n odd. 


Because  the  sign  of  L(y)  is  always  the  same  as  that  of  y\  it  follows  that  the  Bn(y)  are 
expressed  as  a sum  of  terms  which  are  all  positive  or  all  negative.  Now  the  power 
series  for  i^,(y)  and  i^  + ^(y)  converges  rapidly  when  20  > N > |y|  and  moreover  con- 
tains terms  of  only  one  sign.  Furthermore,  the  recursion  relation. 


i 

n - 


,(x)  = (^j-4  in(x)  + in+1(x) 
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when  used  downward  in  n,  does  not  contain  cancellations  so  that  significant  figures 
arc  maintained.  Thus  the  Bn(y)  arc  determined  by  numerically  accurate  methods  for 
all  |y|  N.  For  |y|  > N,  it  was  found  empirically  that  the  closed  form  is  adequate. 

Finally,  the  round-off  error  propagation  was  minimized  by  appropriate  use  of  buffer 
(i.  e.  higher  accuracy)  storage  registers.  Spot  checks  of  the  generated  functions  for 
various  x and  y indicated  6 to  7 significant  figure  accuracy  for  N $ 20. 

Inasmuch  as  it  was  a relatively  simple  modification  of  the  preceding  subrou- 
tine, a second  subroutine  has  been  developed  for  generating  the  set  of  in(x)  where 
n:  0(l)N.  These  modified  spherical  Bessel  functions  occur  in  the  expansions  often 
used  in  the  evaluation  of  multicenter  integrals  and  moreover  inK.  Ru  idenberg's  method 
of  evaluating  two-center  exchange  integrals  of  Slater  AO's.  ^ A general  program  for 
the  latter  method  is  currently  being  coded  for  Whirlwind  by  P.  Merryman  and  K. 
Ruedenberg  as  a part  of  the  general  two-center  integral  evaluation  program  of  Mulli- 
Ken's  group  at  the  University  of  Chicago. 

One-Electron,  Two-Center  Integrals^ 

Since  a large  number  of  one-eiectron,  two-center  integral  values  were  de- 
sired by  the  writer  and  because  such  integrals  occur  often  in  other  molecular  calcu- 
lations, a general  production-style  program  was  written.  This  program  depends 
primarily  on  the  Cn  and  Bn  generation  subroutine  described  earlier  and  allows  the 
calculation  of  one-electron  integrals  involving  wave  functions  which  can  be  expressed 
a linear  combination  of  an  arbitrary  number  of  Is,  2s,  and  2p  Slater  AO's  (e.  g.  a 
fitted  HF  wave  function).  Furthermore,  by  expressing  potential  functions  in  the  form, 


Z (r) 
P 


1^ 

r 


d .r 

l 


the  program  can  also  be  used  to  evaluate  all  of  ihe  potential -type  one-electron  inte- 
grals. The  program  in  its  final  form  can  be  manually  modified  in  a trivial  way  by 
the  computer  operators  to  give  any  of  the  possible  integral  types  (i.  e.  overlap,  kinetic 
energy,  etc.  ).  The  user  gives  as  input  data  a listing  of  the  composition  of  the  wave 
functions  (and  potential,  if  required)  in  terms  of  the  Slater  AO  parameters  and  coeffi- 
cients. The  results  consist  of  a printed  list  containing  the  integral  type,  the  integral 
value,  and  the  wave  function  and  potential  compositions,  so  that  identification  is  com- 
plete and  data  input  mishaps  are  obvious.  A major  advantage  of  the  program  is  that 
the  lengthy  intermediate  results  are  processed  without  error  within  the  machine  and 
need  not  be  brought  out  of  the  computer. 
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Fitting  of  HF  Wave  Functions^ 

Another  application  of  Whirlwind  has  been  a partial  mechanization  of  th>_ 

process  of  fitting  a HF  wave  function  with  exponentials.  The  technique  developed  is 

an  iterative  one  and  is  only  suitable  for  a high-speed  computer.  After  removing 

radial  factors  or  r and  node  factors,  r r,  from  the  tabulated  HF  wave  function, 

o 

there  remains  the  exponential -like  function,  E(r^).  For  an  n-expor.ential  fit,  n pairs 

of  r. 1 and  r-"  are  chosen  from  the  different  regions  of  r^  such  that  r ^ 1 > r^"  > r^1  > 

r,"  > . . . > r 1 > r ".  If  we  define 
2 n n 


Ej(r)  = Z exp(-  b.^  r) 
i 

to  be  the  j1*1  approximation  to  E(r),  then  the  following  prescription  of  cyclical  equa- 
tions can  be  given  where  an  iteration  is  comprised  of  i successively  assuming  the 
values  from  1 to  n.  Defining 


D^(r  /) 


E(r.)  - (1  - XJ) 


Z aj.+  1 exp(-bJlf+  1 v.)  + Z awJ  expf-bjr.) 


k < i 


'k  ‘i'  ’ ^ "k 

k > i 


then 


b-?  = (r . " - r . ')  in 

l i i 7 


rDJ(r.’) 

DJ(r.M)J 


a^  + * = D^(r.')  exp(b^  + 1 r. ') 

l 1 r 1 l 


where 


a ° = b.°  = 0,  0 < X < 1 


Examination  of  these  equations  reveals  that  as  j increases  to  a large  value,  J,  the 
coupling  factor,  1 - X'',  approaches  unity,  and  the  cyclical  equations  express  the  2n 
conditions  of 


EJ(r.')  = E(r.')  and  EJ(ri'')  = EfrT) 

Thus,  if  X is  not  too  much  smaller  than  one  and  J is  large  enough,  the  set  of  a^  and 

b.^  will  smoothly  approach  the  desired  parametric  values.  For  convenience,  after 

1 J 

convergence  is  attained,  the  error  function,  E (r^)  - E(r^),  is  automatically  plotted 

and  photographed  on  the  oscilloscope  camera.  From  the  shape  of  the  error  curve, 

new  choices  of  the  fitting  points,  r.'  and  r.",  can  be  made  and  with  a few  re-runs, 
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optimization  of  the  fit  can  be  essentially  attained.  Empirically  it  appears  that  all  HF 
wave  functions  can  he  fitted  in  this  manner. 

Generation  of  Matrix  Elements 

In  view  of  the  previous  remarks  concerning  function  generation,  it  is  obvious 
that  similar  conclusions  can  be  drawn  concerning  the  generation  of  matrix  elements 
in  secular  equation  problems.  In  particular,  when  using  the  tight  binding  scheme  in 
crystal  problems,  it  is  usually  possible  to  cast  the  matrix  elements  in  the  general 
form  of 


H.  = £ [c^RjKcos  k ■ R,)  + LL.tity  sin  (k  R,)] 

In  writing  a program  for  a specific  crystal  structure,  one  includes  the  properties  of 
the  Cjj(Rj)  and  D.^(R^).  Consequently,  one  need  only  specify  k,  the  Brillouin  zone 
vector,  to  completely  generate  the  desired  matrix  to  be  solved  in  a secular  equation. 
This  procedure  is  found  to  be  especially  convenient  in  the  writer's  current  work  on 
graphite  since  the  S and  H inatr  jes  arising  in  the  non-X  type  secular  equation,  HX  = 
SXX,  both  have  the  same  generation  properties 
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8.  LIMITED  CONFIGURATION  INTERACTION  TREATMENT  OF  THE  NH3  MOLECULE 

The  determination  of  molecular  orbitals  for  the  ammonia  molecule  as  best 
linear  combinations  of  atomic  orbitals  has  been  completed  for  both  the  equilibrium 
case  and  the  planar  case.  Meckler's  mechanization' * ^ of  Roothaan's  formalism^ 
was  used  for  this  process.  The  criterion  of  self-consistency  used  was 


E 

n 


- E 


E 

n 


< 10 


-4 


where  En  is  the  total  molecular  energy  after  the  n**1  cycle.  The  digital  computer 
(Whirlwind)  went  through  approximately  forty  cycles  in  order  to  achieve  consistency. 

The  results  are  generally  quite  gratifying.  For  the  equilibrium  case  the 
calculated  binding  energy,  . 840  a.  u.  , is  about  90  percent  of  the  observed  value,  . 930 
a.  u.  The  hump  in  the  curve  of  molecular  energy  versus  distance  of  the  nitrogen  atom 
from  the  plane  of  the  hydrogens  however  is  wrong  by  an  order  of  magnitude:  calculated 
. 358  a.  u.  , observed  .019  a.  u.  Some  error  in  this  direction  was  expected  due  to  the 
fact  that  the  equilibrium  N - H distance  was  used  for  the  planar  calculation.  However, 
such  a large  error  was  not  foreseen. 

The  one-electron  functions  exhibit  the  expected  behavior.  For  each  symmetry 
type  there  is  one  orbital  somewhat  concentrated  between  the  nitrogen  and  the  hydrogen 
plane  which  is  occupied,  and  one  exhibiting  a node  in  this  region  which  is  unoccupied. 
For  Aj  symmetry  there  is  also  a function  which  is  primarily  directed  away  from  the 
binding  region  (Lennard -Jones  lone  pair)  and  a function  consisting  almost  entirely  of 
nitrogen  1 - s. 

The  results  of  the  calculation  are  now  being  analyzed  in  more  detail.  In  addi- 
tion, preparation  is  being  made  to  use  the  molecular  orbitals  as  the  basis  of  a limited 

(3) 

configuration  interaction  described  earlier. 


1.  A. 

2.  C. 

3.  H. 
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9.  VALENCE  BOND  THEORY.  SIMPLE  SATURATED  MOLECULES 


The  revised  valence  bond  theory,  developed  elsewhere^ ’ ^ and  already  ap- 

(3) 

plied  to  the  iT-electron  systems  of  simple  conjugated  molecules,  is  now  being  ap- 
plied in  uoii-empirical  calculations  on  simple  saturated  molecules.  A preliminary 

treatment  of  the  water  molecule  has  been  completed  and  calculations  on  ammonia  are 

(4  5) 

in  progress.  Both  these  molecules  have  been  studied  previously,  ' using  the  ap- 
proximate self-consistent  field  (SCF)  scheme,  reasonable  binding  energies  being 
obtained.  The  present  calculations  have  a different  object;  they  are  concerned  with 
the  following  questions;  (l)  Can  the  valence  bond  approach  lead  more  directly  to  re- 
sults comparable  with,  or  better  than,  those  of  SCF  MO  theory?  (2)  Can  the  directed 
"hybrid"  orbitals,  commonly  employed  in  empirical  valence  theory,  be  used  to  ad- 
vantage for  this  purpose?  (3)  Is  it  possible  to  use  a model  in  which  some  electrons 
occupy  a "closed  shell"  (composed  of  atomic  inner  shells  and  "lone  pairs")  and  simply 
provide  an  effective  field  for  those  held  responsible  for  binding? 

The  procedure  being  followed  will  oe  illustrated  by  the  case  of  H^O.  The 
guiding  principle  upon  which  the  qualitative  interpretation  of  chemical  valence  and 
stereochemistry  is  usually  based  is  that  of  "maximum  overlapping"  --  crudely  put, 
bonds  are  best  described  by  using  those  orbitals  of  the  different  atoms  which  overlap 
most  strongly,  and  this  principle  is  well  founded  insofar  as  the  use  of  such  orbitals 
provides  a simple  and  direct  way  of  describing  an  accumulation  of  electron  density 
in  the  bond  regions.  In  describing  H^O  then,  where  the  bond  angle  is  considerably 
greater  than  the  90°  between  2p  orbitals,  we  first  imagine  the  oxygen  atom  "promoted’ 
from  the  ground  state  not  to  the  valence  state  (2s)^  (2p  )2  (2p  )*  (2p  ) 1 but  to 

2211  x y z 

u3)  u4)  up  (t2>  where  tj,  are  combinations  of  the  2s  and  2p  orbitals  which 

"point"  towards  the  hydrogen  atoms  Hj  and  (see  Fig.  9-1)  so  as  to  procure  greater 

overlap  with  their  Is  orbitals,  hj  a;  d h2  and  t^,  t^  are  the  res- 
idual orbitals,  determined  by  symmetry  and  orthogonality,  point 
ing  away  from  the  hydrogens  and  describing  the  lone  pairs.  We 
employ  these  orbitals  in  a preliminary  treatment  not  because 
they  are  necessarily  the  most  suitable,  but  because  they  are 
easily  set  up  and  chemically  plausible;  their  deficiencies  can 
always  be  corrected,  if  need  be,  by  configuration  interaction. 
Denoting  oxygen  orbitals  by  S,  S,  X,  Y,  Z,  the  following  symmetry  orbitals 
may  be  defined 

<S,  S,  Z.  o-j(=2'  1/,2(hj  + h2))  ; Y,  <r2(=  Z~  ^(hj  - h2»  : X 

where  o-j,  <r 2 are  left  unnormalized  since  it  is  convenient  to  have  them  related  to  h. 
and  h2  by  a unitary  transformation.  This  basis  of  orbitals,  which  we  shall  derate  by 
the  row  matrix  (A)  = (S,S,  Z,  cr  ^ , Y,  tr2>  X),  is  the  one  used  by  Ellison  and  Shull 
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who  computed  all  the  necessary  integrals.  We  must  now  transform  to  the  basis  (a)  = 

(S  , tj,  t^,  t^,  t^,  hj,  h^):  writing  (A)  = (a)(T),  the  unitary  matrix  (T)  is  very  easily 
determined.  The  straightforward  VB  procedure  would  now  he  to  set  up  a spin  eigen- 
function with  orbitals  t j , and  t^,  h.,  singly  occupied  and  having  "paired"  spins  -- 

this  single  "structure"  representing  the  approximation  of  "perfect  pairing":  this  eourse 

(3) 

is  complicated  by  the  laek  of  orthogonality  of  the  orbitals;  but,  as  we  have  seen,  ' is 
certainly  included  in  a "VB"  treatment  where  the  orthogonalized  orbitals  t t^  --  are 
employed  and  "polar"  structures  (in  which  orbitals  may  be  empty  or  doubly  occupied) 
are  admitted  and  the  VB  calculation  can  be  carried  out  without  the  usual  approximations 
of  VB  theory.  We  therefore  form,  from  (a),  a set  of  orthonormalized  orbitals  (a),  by 
the  Lowdin  method,  i.  e.  , (a)  = (a)(S  ) Here  (S  ) is  the  matrix  of  overlap  inte- 

grals  (metrical  matrix)  for  the  (a)  basis  and  is  related  to  that  for  the  (A)  basis  by  (S^)  = 
(T)(S  . )(T)  since  (T)  is  unitary  (S  ) = (T)(S.)  */^(T)^  and  it  follows  that  the 

orthonormalized  AO's  (directed  AO's  included)  are  related  to  the  symmetry  orbitals  by 
(a)  = (A)(p)  where  (p)  = (S^)  '^(T)1.  All  integrals  referred  to  the  (a)  basis  may  now 

be  expressed  in  terms  of  those  given  by  Ellison  and  Shull. 

The  next  step  in  the  VB  calculation  is  to  set  up  suitable  VB  structures  in 
terms  of  the  AO's.  First  we  limit  ourselves  to  those  in  which  S,  t^  and  t^  are  doubly 
occupied,  and  we  deal  explicitly  only  with  the  electrons  occuping  t , t^  Ej  and  h^,  us- 
ing an  effective  one-electron  Hamiltonian,  h,  which  takes  into  account  the  closed  shell, 
instead  of  that  for  an  electron  in  the  field  of  the  bare  nuclei  (f,  say).  This  is  a quite 
legitimate  reduction  (see  Ref.  2)  so  long  as  we  keep  to  configurations  showing  a filled 
elosed  shell  and  leads  effectively  to  a four-electron  problem.  The  effective  Hamilton- 
ian h may  be  defined  by  its  matrix  elements  between  two  of  the  valenee  orbitals,  p and 
q say: 


(p|h|q)  - (p  | f | q)  + 2£(xp|g|xq)  - £ (xp  | g | qx) 

x x 

where  x is  a closed  shell  orbital  and  summations  are  over  all  orbitals  of  the  elosed 
shell  (,S  i tj.  t^).  The  energy  calculated  from  the  four-electron  problem  is  then  relative 
to  a elosed  shell  or  "framework"  energy 

Ef  = 2 Y (x  | f | x)  + Y.  (xx  |g|  xx)  + 4 Y (xy|g|xy)  - 2 £ (*yjg|yx) 
xx  xy  xy 

summations  again  being  over  the  closed  shell  (eour.ting  distinct  pairs  xy).  The  four- 
electron  VB  structures  employed  so  far  are  shown,  with  an  obvious  symbolism,  in 
Fig.  9-2.  The  first  five  are  the  "obvious"  ones  which  would  arise  as  the  dominant 
terms  in  expansion  of  the  perfeet  pairing  approximation  into  VB  structures,  and  which, 
in  the  eharge-hopping  interpretation  should  be  capable  of  describing  the  two  O-H  bonds. 


A 


-32- 


aft* 


(VALENCE  BOND  THEORY.  SIMPLE  SATURATED  MOLECULES) 


/ \ ./  \ / V \ / \+ 
_L.  -/  V +/  V 


The  remainder  correspond  to  intra-atomic  atomic  spin  pairings  and  double  "excitations" 
but  either  have  low  "energies"  (diagonal  matrix  elements)  or  interact  strongly  with  the 
first  five.  Thus,  the  "resonance"  integral  (p|h|q),  is  associated  with  binding  in  the 
overlap  region  of  p and  q containing  the  interaction  energy  of  the  nuclei  and  a charge 
distribution  pq  in  this  region  and  governs  the  interaction  between  structures  which  show 
a charge  shift  between  p and  q.  As  anticipated,  it  is  very  large  for  the  strongly  over- 
lapping orbitals  of  the  bends,  but  is  small  otherwise  (e.  g.  when  p and  q belong  to  differ- 
ent bonds),  and  consequently  structures  showing  intra-atomic  charge  shifts  (Fig.  9-2) 
or  crossed  links  should  participate  only  weakly,  doing  little  to  stabilize  the  molecule. 

But  the  doubly  polar  structures  are  likely  to  be  important  --  though  we  are  still  omitting 
those  showing  a "doubly  ionized"  oxygen  atom,  which  lie  rather  higher  in  energy. 

This  general  picture  is  supported  by  the  actual  calculations.  Unfortunately, 
detailed  numerical  results  are  not  on  hand  (the  writer  being  at  present  en  route  to  King's 
College,  Newcastle-on-Tyne,  England)  but  the  following  qualitative  conclusions  are 
worthy  of  note. 

There  is  no  difficulty  in  handling  the  VB  calculation  by  the  methods  already 
developed.  The  only  real  computational  obstacle  is  the  transformation  to  integrals 
over  the  AO's,  and  this  is  easily  mechanized^  being  performed  on  Whirlwind  in  a 
very  few  minutes  once  the  (28  X 28)  matrices  have  been  written  down.  * The  iterative 
SCF  procedure  is  replaced  by  the  solution  of  a final  secular  determinant  --  again  easily 
mechanized  --  and  there  is  little  doubt  that  the  accuracy  of  the  SCF  result  could  be 
surpassed  without  too  much  trouble  by  adding  a few  more  structures  (e.  g.  breaking 
into  the  closed  shell).  But  the  whole  object  of  these  calculations  was  to  find  a simple 
way  of  setting  up  good  wave  functions,  making  use  of  the  basic  chemical  concepts  of 
directed  valence  and  lone  pairs  and  without  any  appeal  to  existing  knowledge  (e.  g.  the 
SCF  results  themselves),  and  from  this  point  of  view  the  results  are  not  altogether 
encouraging.  The  decrease  in  total  electronic  energy,  accompanying  molecule  forma- 
tion, is  found  to  be  about  95  percent  of  that  predicted  using  the  full  SCF  method  (and 
the  wave  function  might  therefore  be  reasonably  goodh  but  this  decrease  is  completely 
offset  by  the  increase  in  internuclear  repulsion  energy  and  the  last  5 percent  is  found 

I am  indebted  to  Dr.  Meckler  for  the  use  of  his  machine  program  for  inis  purpose. 
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to  be  crucial  in  getting  a binding  ener  gy.  This  is  a familiar  situation  in  binding  energy 
calculations  --  the  electronic  energy  is  good  but  the  binding  energy  is  very  poor.  The 
reason  for  this  failure,  however,  does  not  lie  wholly  in  the  inadequacy  of  the  VB  struc- 
tures selected,  for  at  the  outset  we  msde  a somewhat  arbitrary  choice  of  orbitals, 
upon  which  to  base  the  limited  configuration  interaction  calculation.  Indeed,  it  is  likely 
that  the  use  of  differently  directed  orbitals  (describing  "bent"  or  "strained"  bonds) 
would  have  led  to  an  improvement  of  the  required  order  still  within  the  framework  of 
a closed-shell  model  with  a simple  chem'oal  interpretation  (inner  shell,  two  lone  pairs, 
two  bonds).  This  belief,  that  the  basic  orbitals  themselves  are  considerably  at  fault, 
was  strengthened  by  making  a limited  SCF  MO  calculation,  using  one  determinant  based 
on  the  configuration  (5)^,  (t^)*,  (t^)^,  (T^)^,  (T^)^,  where  T f = N(tj  + \hj)  and  T ^ = 

N (T2  + \h2). 

The  only  parameter  admitted  --  X,  which  describes  the  polarity  of  the  O-H 
bond  and  upon  which  the  energy  is  quite  strongly  dependent  --  was  directly  varied  and 
the  final  minimum  value  was  found  to  be  no  better  than  that  reached  in  the  VB  calcula- 
tion. 

These  results  therefore  suggest  thai  a simple  approach,  in  which  attention  is 
focussed  upon  the  individual  bonds,  cannot  be  quantitatively  successful  unless  the  bond 
orbitals  (localized  MO's)  of  an  MO  treatment  or  the  directed  orbitals  of  a VB  (or  VB) 
treatment  are  chosen  with  extreme  care.  Unless  this  can  be  done  (as,  in  effect,  it  is 
in  a full  SCF  treatment  where  none  of  the  orbitals  are  chosen  arbitrarily)  there  seems 
to  be  no  alternative  but  to  employ  a "hammer  and  tongs"  approach  in  which  more  nd 
more  configurations  are  thrown  in  to  correct  for  an  inappropriate  choice.  There  can 
be  no  doubt  that  this  procedure  is  feasible  (many  more  configurations  were  handled  in 
the  calculations  on  benzene^)  and  the  final  results  can  always  be  given  a simple  pic- 
torial interpretation  by  shifting  attention  onto  the  actual  charge  density  (which  can  be 
calculated  in  many-configuration  theories  by  the  methods  described  elsewhere!'*’  ^ 

But  this  is  a clumsy  way  of  getting  a solution  and  it  is  to  be  hoped  this  procedure  car. 
be  largely  obviated  by  using  intuition  and  previous  experience  in  choosing  more  satis- 
factory basic  orbitals.  Again,  it  seems  that  a semi-empirical  development  is  suitable; 
for  binding  is  primarily  associated  with  a few  large  integrals  and  these  might  well  be 
regarded  as  disposable  constants  and  the  basic  orbitals  themselves  would  then  no 
longer  play  a crucial  part  in  the  calculations. 

The  writer  would  like  to  take  this  opportunity  of  recording  his  gratitude  to 
Professor  Slater,  and  to  members  of  the  Solid-State  and  Molecular  Theory  Group,  for 
the  hospitality  he  has  enjoyed  at  the  Institute  and  for  the  benefit  of  many  discussions 
during  a most  happy  and  stimulating  year. 
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10.  NATURAL  ORBITALS  FOR  THE  HELIUM  ATOM 


In  his  contribution  to  the  last  Quarterly  Progress  Report,  P.  O.  Lowdin^ 
introduces  the  concept  of  natural  orbitals  and  discusses  some  of  their  interesting 
properties.  He  introduces  them  in  the  following  manner.  We  imagine  that  we  are 
trying  to  approximate  the  solution  of  an  n electron  problem  with  a Hamiltonian  con- 
taining the  usual  kinetic  and  potential  energy  terms.  In  order  to  approximate  the  solu- 
tion to  this  problem,  we  form  antisymmetric  stat  s formed  by  taking  linear  combina- 
tions of  Slater  determinants  each  containing  n one-electron  spin  orbitals  selected 
from  N spin  orbitals,  j . . . 4^.  If  we  call  the  wave  function  described  above  J (xj  .. 

. . x ),  a first  order  density  matrix  can  be  defined  as 


v(x  j ' lx  | ) = N J ^(x,1,  x2  . . . xn)  $(xr  x. 


x ) dx , 
n'  2 


dx 


(10. 


This  density  matrix  can  be  written  in  terms  of  our  one-elcctror.  spin  orbitals  in  the 
form 

, N 

yUj'Ixj)  = Z (k.i)  \ (xj  ')  tp/x  1 ) Y^  (10.2) 

1 


Here  tile  constants  y^  consist  of  some  forms  in  the  constants  that  multiply  the  de- 
terminants in  the  wave  function  J.  If  we  go  to  another  basis  in  the  one-electron  func- 
tions by  a unitary  transformation,  we  shall  induce  a unitary  transformation  on  the  co- 
efficients yi;  • Let  us  call  the  transformed  functions 
•ik 

<t>k  = Z(«)  +Q  oak  (10-3) 

W'e  can  now  look  for  the  transformed  density  matrix  which  is  diagonal.  That  is  we  can 
choose  our  unitary  transformation  such  that 

U^yU  ” n = diagonal  matrix 

yUj'Uj)  = Z(R)  nk  \*(xl  ')  \(xi)  (10.4) 

We  can  see  that  in  general  the  natural  spin  orbitals  will  mix  the  original  spin 
orbitals  of  a and  p spin.  Under  the  simplifying  assumption  that  our  total  wave  func- 
tion ^ has  Mg  (total  z component  of  the  spin)  as  a good  quantum  number  this  will  not 
be  the  case  since  the  density  matrix  y^  v.ill  be  blocked  off  into  two  parts.  The  one 
part  we  might  denote  as  y*k  and  the  other  part  we  shall  call  y^.  The  constants  y^ 
are  the  coefficients  of  the  products  of  orbitals  with  a spin  in  (3.  2)  and  the  y^,  are 
the  coefficients  of  products  of  orbitals  with  p spin.  There  are  no  cross  terms  in  this 
case.  The  reason  for  this  can  be  seen  by  considering  the  origin  of  any  term  in  the 
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expansion  of  the  density  matrix  ( 1 0.  2).  Any  term  in  the  sum  in  (10.  2)  will  arise  from 
the  integration  of  the  product  of  two  Slater  determinants  over  n - 1 coordinates.  If 
^k(xi')  has  Q spin  and  has  p spin  the  integration  over  the  coordinates  x^  . . . >-n 

will  yield  zero  since  the  integral  over  spins  will  be  the  integral  of  the  product  of  two 
functions,  one  having  spin  Mg  + 1 / 2 and  the  other  having  total  spin  Mg  - 1/2.  Since 
our  density  matrix  is  blocked  off  into  two  parts,  one  for  a and  the  other  for  p spin,  the 
unitary  transformation  which  diagonaliz.es  this  matrix  will  also  be  blocked  off.  We 
shall  thus  have  no  mixing  of  a and  p spin  orbitals. 

Part  of  the  hope  that  Ldwdin  expressed  in  connection  with  the  natural  orbitals 
was  that  a configuration  interaction  treatment  of  the  problem  of  the  many-eiectron  sys- 
tem would  be  more  rapidly  convergent  in  terms  of  the  natural  orbitals.  In  order  to 
study  the  natural  orbitals  in  a simple  case  and  to  study  the  convergence  of  the  con- 
figuration interaction  in  terms  of  natural  orbitals,  it  was  decided  to  carry  out  a simple 
calculation  on  the  helium  atom. 

For  one-electron  orbitals  a Is  and  a 2s  wave  function  were  chosen  in  the  form 

u,  . [z/ y47u')3/2J  e- Z'r  ^ 

u?  = [l/uyryuo]  (z*)3/2  (2  - Z')  e 2 (10.5) 

By  making  all  possible  *S  states  for  this  two-electron  problem  we  arrive  at  three  pos- 
sible configurations:  (ls)Z,  (ls2s),  (2s)^.  Our  many -electron  wave  function  is  now 
expressed  as  a sum  of  three  terms 

J(xj,x^)  = c^ls)2  + c.,(ls2s)  + c3(2s)2  (10.6) 

The  expectation  value  of  our  Hamiltonian  was  taken  with  respect  to  this  wave  function 
for  a given  value  of  the  effective  charge  Z'  and  the  values  of  the  energy  and  the  c's 
we  •«  determined  ^ It  was  found  that  the  best  value  of  Z'  was  roughly  28/  16  and  that 
for  this  value  the  c's  and  the  energy  were 

E = - 5. 70939 
Cj  = .993173 

cz  = - . 115706 

c,  = - . 014819  (10. 7) 

Table  10-1  contains  a comparison  of  the  energy  obtained  in  this  way  with  the  Ilartree- 
Fock  energy,  the  experimental  energy  and  the  best  possible  energy  obtainable  from  a 
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(is)  configuration  with  the  wave  function  of  the 


type  (10.  5). 


Table  10-1 

Energy  of  the  helium  atom  by  various  approximations 
(energies  in  Rydbergs) 


Experimental 
llartrce  -Fock 
(is)2  Z'  = 27/16 

(Is)2  + ( 1 s2s)  + (2s)2  Z'  = 28/16 


- 5. 80736 

- 5.  75 

- 5. 6953 

- 5. 709394 


Knowing  the  wave  function  obtained  in  this  manner  it  is  possible  to  construct 
a density  matrix  for  this  problem. 

y(x  j|  x j ')  = . 993085  iijd')  a(l')  Ujd)  a(l)  + . 0069  1 4 u^(  1 ')  a(l ')  t>2(  1 ) a(l ) 

- . 080045  Ujd1)  a(l')  u2(l)  a(l)  - . 080045  u2(  1 ')  n(  1 ')  u ,( 1 ) a(l ) 

+ identical  terms  with  all  a's  replaced  by  p's.  (10.  8) 


Therefore 

v|j  = Yn  = -993085 

ylZ  = yZl  = V12  = Y21  = • 080045 

y\z  = y\7  = ■ 006914  ( 10.  9) 

Vie  notice  here  the  factoring  of  the  density  matrix  into  a plus  spin  part  and  a minus 
spin  pari  which  we  mentioned  earlier.  We  also  aotice  that  the  plus  and  minus  spin 
parts  of  the  density  matrix  are  identical.  This  arises  from  the  fact  that  the  state  we 
are  considering  is  a singlet  and  therefore  the  density  matrix  must  remain  unchanged 
if  we  replace  all  a's  by  p's.  In  order  to  find  the  natural  orbitals  all  we  need  da  is  to 
diagonalize  the  two  by  two  matrix 
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We  find  that  if  we  diagonalize  this  matrix  we  get  as  the  diagonal  elements  in  the  trans- 
formed matrix 


r.j  = .999541 

n^  = . 000459  (10. 11) 

The  transformation  which  accomplishes  this  diagonalization  has  elements 


U11  - 

. 996764 

U22  = 

. 996764 

U12  = - 

. 080379 

u,,  = + 

<_  l 

. 080379 

(10. 12) 

Therefore  our  natural  orbitals  are  given  by 

Uj'  = . 996764  u.  - . 080379  u2 

u2'  = . 080^79  + . 996764  u2  (10.13) 


The  new  Is  function  therefore  has  added  to  it  a small  amount  of  2s  and  the  new  2s  func- 
tion has  added  to  it  a small  amount  of  Is. 

We  are  now  in  a position  to  determine  the  convergence  of  the  configuration 
interaction  treatment  in  terms  of  the  natural  and  the  original  orbitals.  We  can  con- 
struct from  our  natural  orbitals  configurations  similar  to  those  which  we  can  construct 
from  our  original  orbitals.  We  shall  call  these  (is1)2  (ls'2s')  and  (2s')“.  We  can  now 
do  our  configuration  interaction  in  three  approximations.  We  can  compare  the  ener- 
gies of  single  (is)  configurations  made  out  of  natural  and  our  original  orbitals,  we 
can  compare  the  energies  we  get  :‘ram  solving  the  two  by  two  connecting  the  (is)2  and 
(ls2s)  and  we  can  finally  compare  the  energies  that  we  get  from  the  complete  three  by 
three  done  in  both  manners.  In  Table  10-2,  we  show  the  results  of  tins  comparison. 

Table  10-2 

Comparison  of  the  energies  of  the  helium  atom  arrived  at 
using  natural  orbitals  and  the  original  orbitals 
(energies  in  Rydbergs) 


(is)2 

- 5. 707462 

- 5. 

6875 

(Is)2  + 

( ls2s) 

- 5. 707531 

- 5. 

708464 

(Is)2  + 

(ls2s)  + (2s)2 

- 5. 709394 

- 5. 

709394 
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We  notice  from  this  table  that  the  natural  orbitals  give  a much  better  energy  for  a 
single  configuration  than  do  the  original  orbitals.  Adding  the  next  configuration  gives 
a very  small  improvement  in  the  energy  in  the  natural  orbital  framework.  For  two 
configurations,  however,  it  would  have  been  advantageous  to  use  the  original  orbitals. 
For  three  configurations  it  makes  no  difference  which  one-electron  functions  we  use 
as  the  basis.  From  this  simple  example  it  is  not  clear  whether  we  have  improved  the 
convergence  by  using  the  natural  orbitals.  For  one  configuration  it  is  better  to  use 
the  natural  orbitals,  for  two  configurations  it  is  better  to  use  the  original  orbitals, 
for  three  i*.  makes  no  difference  which  we  use  since  we  ar  e doing  the  complete  con- 
figuration interaction  which  can  be  carried  out  using  our  basic  set  of  one-electron 
orbitals. 

Even  the  result  that  the  one  configuration  ener  gy  is  improved  by  the  use  of 
natural  orbitals  is  an  uncertain  result  as  we  can  show  by  a simple  argument.  Imagine 
that  ne  did  our  configuration  interaction  problem  for  helium  using  the  Hartree-Fock 
Is  and  2s  functions  as  our  basic  set  of  functions  rather  than  the  ones  in  (3.  5).  In  *his 
case  we  would  certainly  have  a three  by  three  configuration  interaction  problem  and 
could  once  again  find  the  natural  orbitals.  In  this  case  the  natural  orbital  would  turn 
out  to  be  some  mixture  of  the  Is  and  2s  Hartree-Fock  functions.  If  we  now  used  the 
natural  orbitals  formed  in  this  manner  to  form  the  (is)2  configurations,  we  could  ask 
if  the  energy  of  the  singlet  state  would  be  improved  by  the  use  of  the  natural  orbitals. 
It  would  certainly  not  be  improved  since  the  original  Hartree-Fock  functions  were 
defined  in  such  a way  that  they  gave  the  best  possible  energy  which  can  arise  from  a 
single  determinant.  An  admixture  of  2s  to  the  Hartree-Fock  Is  could  only  have  the 
effect  of  raising  the  energy. 

From  the  simple  example  and  the  discussions  it  would  seem  that  it  is  doubt- 
ful that  natural  orbitals  are  going  to  increase  the  speed  of  convergence  of  a configura- 
tion interaction  approach  to  the  many-electron  problem. 
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11.  HARTREE-FOCK -SLATER  SELF-CONSISTENT  FIELD  FOR  Fe  (3d)6  (4s)2 

Solution  of  the  one -electron  wave  functions  for  Fe  with  an  unbalanced  spin 
configuration,  namely  five  3d  electrons  with  spin  up,  one  3d  electron  with  spin  down, 
and  two  4s  electrons  (one  of  each  spin)  has  reached  self-consistency  according  to  the 
criterion  established  by  Hartrec.  The  3d  and  4s  functions  show  the  most  striking 
variations.  In  both  cases,  the  spin-up  function  is  shifted,  with  respect  to  the  spin- 
down  function,  toward  the  origin  (see  Table  11-1).  Both  the  3d  and  4:;  electrons  spend 
considerable  time  in  regions  where  the  averaged  exchange  potential  is  a sizable  frac- 
tion of  the  total  potential  function  appearing  in  the  radial  wave  equation.  Thus,  the 
two  exchange  potentials  become  important  in  determining  the  behavior  of  the  radial 
wave  functions  for  3d  and  4s. 


Table  11-1 


p 

r 

r(M . G.  ) 

1 s+ 

. 035 

. 037 

ls- 

. 035 

2s+ 

. 225 

2s- 

. 225 

. 221 

2p+ 

. 180 

2p- 

. 180 

. 181 

3s+ 

. 700 

3s- 

. 705 

. 735 

3p+ 

. 705 

. 735 

3p- 

.715 

3d+ 

. 680 

. 735 

3d- 

. 700 

4s+ 

2.  240 

4s- 

2.  360 

2.  697 

Tile  charge  density  for  the  other  functions  (Is,  2s,  2p,  3s,  3p)  is  essentially 
the  same  with  respect  to  spin  interchange.  The  above-mentioned  shift  is  not  prominent, 
althuugh  still  present.  The  total  potential  seen  by  these  electrons  is  largely  made  up 
of  ordinary  plus  the  angular  momentum  term,  the  exchange  potentials  being  here 
less  important. 

In  Table  11-2  we  list  the  radii  of  the  electronic  orbits  (maxima  of  p2(r))  for 
our  final  self-consistent  functions  and  for  those  determined  by  Manning  and  Goldberg. 
This  table  is  of  limited  value,  but  it  serves  to  indicate  the  general  trend  of  our  calcu- 
lations. The  exchange  has  tended  to  pull  the  charge  distribution  in  toward  the  origin 
and  to  lower  the  one -electron  energies  as  one  might  guess.  However,  starting  with 
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Table  11-2 

Energies  (- Rydbergs) 


M 

.G. 

Present  Calculations 

Ionization  Potentials 

! s 

523.  0 

534.  5 (i) 

584. 2 (-) 

523.  9 

2 s 

60.  79 

61.  19  (+) 

60.  67  (-) 

62.  5 

2p 

53.  02 

57.  09  ( + ) 

55.  90  (-) 

52.  2 

3s 

6.  973 

7.  463  (+) 

6.  930  (-) 

6.  9 

3p 

4.  600 

5.  061  (+) 

4.  540  (-) 

4.  1 

3d 

0.  7578 

1 . 1223  (-1) 

0.  6636  (-) 

0.  60 

4s 

0.  4836 

0.  5315  ( + ) 

0.  4275  (-) 

0.  58 

3s  in  the  energy  table,  the  spin-down  functions  have  a higher  energy  than  do  the 
Manning-Goldberg  functions. 

The  inward  shift  of  the  one-electron  radial  charge  densities  is  more  marked 
for  the  functions  starting  with  3s,  and  the  amount  of  shift  seems  to  follow  an  argu- 
ment like  that  given  for  the  differences  between  the  spin-up  and  spin-down  functions; 
one  is  on  less  firm  ground  here,  however,  in  view  of  the  self-consistency  of  the  two 
sets  of  calculations  (Manning-Goldberg  and  present  calculations). 

J.  H.  Wood 
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12.  NUCLEAR  ELECTRIC  QUADRUPOLE  INTERACTION  IN  THE  KC1  MOLECULE 


The  calculation  of  polarizability  effects  in  F is  being  carried  out  in  the  re- 
vised scheme  discussed  in  a previous  Progress  Report.  ^ Unfortunately,  this  new 
formulation  has  rearranged  the  problem  to  such  an  extent  that  it  has  been  necessary 
to  re-evaluate  the  integrals.  However,  the  experience  gained  previously  has  been 
helpful  and  work  is  proceeding  rapidly. 

Refe  pence 
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13.  ELECTRON-LATTICE  INTERACTIONS 


The  problem  of  interaction  between  conduction  electrons,  or  impurity-bouod 
electrons  and  lattice  vibrations  where  the  interactions  arc  too  strong  to  be  handled  by 
perturbation  methods  (e.  g.  ionic  crystals)  has  been  approached  in  a variety  of  differ- 
ent ways.  Frohlich,  Pelzer,  and  Zienau*  and  Gross^  have  worked  in  a Fock  space 
for  the  phonon  assembly  inti  ouucing  cutoffs  at  one,  two,  or  three  quanta,  then  attempt - 


(3) 


mg  an  exact  solution  of  'he  ?nnrnvimate  problem  remaininn.  Lee,  Low,  and  Pines, 

(4) 

and  Gurari  ' have  applied  intermediate  coupling  theory  in  what  amounts  to  a Hartree- 
Fock  approach  in  the  momentum  representation  of  Fock  space.  Various  Russian 
writers  including  Bogolyubov,  ^Tyablikov,  and  Pekar^^  have  worked  in  an  electron- 
lattice  configuration  space  making  various  approaches  to  an  adiabatic  approximation. 

Basic  to  all  these  approaches  are  certain  assumptions.  The  periodic  poten- 
tial of  the  lattice  is  taken  into  account  through  the  introduction  of  an  effective  mass. 

The  lattice  is  treated  as  a harmonically  vibrating  continuum  of  polarization  charge 

(8) 

described  by  a phenomenoiogica:  Hamiltonian  introduced  by  Frohlich.  ' The  analyses 
are  restricted  to  small  momenta,  that  is  slow  electrons. 

(9) 

Several  quest. ons  have  been  raised:  Haken  in  treating  a highly  specialized 
model  has  suggested  that  the  effect  of  the  actual  periodic  potential  on  the  effect  of  in- 
teraction may  be  considerable.  Lee  and  Pines  claim  to  obtain  an  exact  solution  in  the 
limit  of  strong  coupling,  yet  Pekar  using  a variational  method  seems  to  obtain  lower 
energies  for  sufficiently  strong  coupling.  Lastly,  the  Fock-space  methods  suggest  an 
effective  mass  increasing  sharply  with  momentum,  but  these  methods  are  valid  only 
for  small  momenta.  Attempts  are  now  being  made  to  resolve  these  tm'-ertainties. 

First  an  unsuccessful  attempt  has  been  made  to  investigate  higher  momenta 
of  the  electron-lattice  system  by  a variational  calculation  in  the  Fock-space  momentum 
representation  in  which  the  correlation  between  different  phonons  was  explicitly  in- 
cluded in  the  trial  function.  The  Frohlich  Hamiltonian  for  an  electron  described  by 
variables  q,  p and  a system  of  traveling-wave  oscillators,  the  k having  frequency 
creation  operator  ak*  and  annihilation  operator  ak  is 

2 
3 

2m 


H 


ik  • q , * 

ake  +Vak 


= £-  + Z gk  a,,  e"'  H + g_*  a,.*  e 

? m % ^ 


ik  • q 


I* 


“k  ak*  ak 


neglecting  zero-point  energies.  Introducing  Jf?  = U ' HU  by  the  unitary  operator 


U = e 


iti  Z k"  ak*  ak  • q" 


one  obtains 
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(;7  - ti  £ iT  ak*  aR) 


2m 


£ 6k  ak  + Sk*  ak*  + I *“k  ak 
k k 


in  which  q is  absent  and  its  conjugate  p,  now  obviously  the  total  momentum,  is  con- 
served. Restricting  ourselves  to  states  of  definite  total  momentum  p.  wc  may  con- 
sider the  operator  p as  a c -number.  If  4n(k  j • • • k^;  p)=  <C  k j ...  k^  p | U 1 
is  the  probability  amplitude  that  a system  in  the  state  |U  1 $>  has  n quanta  of  mo- 
menta tikj,  Tik^  , . , hk^  respectively  and  total  momentum  p,  then 


<i|H-E|J>=  I I I ■ ■ Z {^n(irr  - - lit)2  +^I  <^k.-  E] 

n kj  k^  kn  Lim  n i i l 

+ v/n  + 1 g(k. ) 4 *(k  . . . . k ; p)  4 . (k . . . . k k;  p ) 

v 6'  r m 1 n ^ ' *n + 1 I n i * ' 

+ Vn  r,  *(k . . . . k ; p ) g*(k  ) 4 . (k  . . . k . ; p )} 

v Jn  1 n H ' 6 n'  *n  - 1 1 n - 1 K ') 


= I(E;  pi 


The  requirement  that  HJ  = EJ  is  equivalent  to  6r  =0  for  arbitrary  variations  sym- 
metric in  all  k As  in  the  infinite  set  of  functions  j • 

The  method  of  FPZ  was  to  set  4 - 0 for  all  n > m,  m = I and  to  solve  the 

(2)  n 

two  remaining  e:;  '?ticns.  Gross  ' attempted  an  extension  to  m = 2 and  3 and  was  led 
to  m + I coupled  integral  equations.  This  method  has  the  advantage  of  including  cor- 
relations between  two  or  three  phonons  accurately,  but  it  is  a good  approximation  only 
if  the  probability  of  rrmre  than  m phonons  being  excited  is  very  small.  Furthermore, 
for  more  than  one  phonon,  the  integral  equations  cannot  be  solved  exactly,  and  for 
more  than  two  even  the  approximate  methods  become  extremely  cumbersome.  Lee, 
Low,  and  Pines  and  Gurari  have  specified  that 


^(k> 


kn;  p)  = 


n f-(k,) 
i=i  p 1 


where  C is  a normalization  constant,  and  Lee  and  Pines  have  generalized  this  ap- 
proach somewhat.  Nevertheless,  the  wave  functions  thus  derived  do  not  explicitly  in- 
clude correlation. 

It  was  suggested  by  llylleraas's  calculations  on  helium  to  set  up  a Fock- 
space  momentum  representative  (40(p  ),  4[(kj:p),  ^(kj,  k-,;  p )...)  in  which  the 
form  of  each  function  4n(kj  • • • k^;  p)  is  specified  but  depends  on  a few  arbitrary 
parameters.  The  functional  I(E;  jT)  becomes  a function  of  these  parameters  with 
respect  to  which  it  is  minimized.  It  was  hoped  that  if  a set  of  functions  could  be 
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found  whose  form  was  suggested  by  Gross's  accurate  wave  functions  for  the  one-  and 

two-quantum  cutoffs  yet  which  would  lead  to  tractable  integrals,  the  minimization 

procedure  could  be  performed  by  the  Whirlwind  digital  computer  for  a cutoff  of  as  high 

as  30  quanta.  Functions  of  the  form 

1 


have  'he  proper  symmetry,  explicitly  include  the  correlation  and  lead  to  multi- 
dimensional integrals  which  can  be  evaluated  analytically,  at  least  for  zero  total  mo- 
mentum. Preliminary  hard  computations  indicated  that  though  some  improvement 
could  be  obtained  if  p.  were  allowed  to  vary  from  zero  (case  of  no  correlation),  never- 
theless the  depressions  of  energy  due  to  the  interaction  were  only  about  half  those 
obtained  by  olher  methods  due  to  the  poor  fit  of  the  one  quantum  trial  function  to  the 
"exact"  one  quantum  function.  Attempts  at  constructing  trial  functions  better  approxi- 
mating the  functions  found  by  Gross  led  to  seemingly  impossible  many-dimensional 
integrals  for  any  significant  number  of  quanta. 

The  fundamental  difficulty  is  that  one  is  trying  to  fit  a series  of  unknown 
functions  rather  than  just  one  as  in  the  many-electron  problem.  If  an  undetermined 
linear  combination  of  m prescribed  functions  is  constructed  to  approximate  each 
then  tiie  matrix  to  be  diagonalized  for  n-quanta  cutoff  is  roughly  of  order  m x n 
so  that  n must  be  quite  small.  If  instead,  more  complicated  single  functions  are 
used  for  a better  approximation,  the  integrals  become  impossible.  The  work  of  Lee 
and  Pines  has  perhaps  gone  as  far  as  possible  by  not  introducing  correlation  into  the 
trial  functions  and  thus  retaining  simple  integrals  regardless  of  the  number  of  quanta. 
Further  attempts  along  this  line  have  been  abandoned,  at  least  temporarily,  since  the 
dependence  of  energy  on  the  correlation  parameter  and  on  the  total  momentum  are 
meaningless  if  the  trial  functions  at  best  are  poor  approximations.  Other  methods  of 
extrapolating  to  high  momenta  are  being  considered. 

A second  investigation  now  under  way  is  stimulated  by  the  desire  to  construct 
a model  which  can  be  solved  with  sufficient  exactness  to  serve  as  a yardstick  in  de- 
termining the  validity  of  the  various  adiabatic  approaches  and  the  validity  of  account- 
ing for  the  periodic  potential  through  an  effective  mass.  Gross*'  * ^ has  determined 
the  exact  solutions  to  the  eigenvalue  problem  arising  when  a free  electron  is  allowed 
to  interact  with  only  one  traveling  lattice  wave  and  all  motion  is  confined  to  one  dimen- 
sion. For  sufficiently  strong  coupling  he  finds  the  state  with  minimum  energy  has 
non-zero  momentum.  Since  the  one-phonon  analysis  is  essentially  spatially  asym- 
metric, the  problem  of  an  electron  with  two  degenerate  modes  (right-  and  left-going) 
in  one  dimension  has  been  set  up  and  programmed  for  Whirlwind  by  Meckler,  using 
a cutoff  of  six  quanta.  An  investigation  of  the  energy-momentum  curves  for  mtermedi- 
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ate  and  $t  rong  coupling,  the  dependence  on  lattice  wave  length  of  the  effect  of  the  in- 
teraction, ana  the  convergence  of  the  cutoff  procedure  are  being  investigated.  It  is 
hoped  to  push  the  one-  or  two-cscillator  model  as  far  as  possible  with  regard  to  the 
per  iodic  potential  and  the  validity  of  the  adiabatic  approximation.  The  Fock-space 
approaches  are  not  applicable  here,  depending  as  they  do  on  the  presence  of  a large 
number  off  independent  oscillators. 
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14.  THERMAL  VIBRATIONS  IN  Cu-Zn  SYSTEM  CRYSTALS 


The  work  described  in  the  previous  report  is  nearing  completion,  and  it  is 
anticipated  that  a paper  will  soon  be  written  on  this  investigation  of  the  atomic  force 
constants  in  copper,  in  addition  to  a Ph.  D thesis.  The  expressions  for  the  change  in 
conduction  electron  charge  density  (Ap)  produced  by  the  displacement  of  an  ion  core 
have  been  extensively  investigated.  Ap  has  been  calculated  at  a number  of  points  in 
the  lattice  as  a function  of  edge-width  of  the  finite  lattice  used,  both  for  cubic  and 
for  narallelopiped  macroscopic  shapes.  The  convergence  shown  is  of  the  nature  of 
damped  oscillation  about  a constant  value  as  asymptote,  and  the  convergence  pattern 
is  the  same  at  different  points,  except  far  from  the  origin,  where  spurious  structure 
appears.  For  two  lattices  extensive  plots  of  Ap  as  a function  of  position  in  the  lattice 
have  also  been  made,  and  found  in  general  agreement.  A lattice  of  ten  atomic  planes' 
width  has  been  settled  upon  for  final  calculations,  as  the  best  compromise  between 
accuracy  and  ease  of  calculation.  The  numerical  Coulomb  integration  of  Ap  to  get  its 
contribution  to  the  atomic  force  constants  is  now  underway,  and  appropriate  tables 
have  been  constructed  for  comparison  of  the  results  with  experimental,  nuclear  dipole, 
closed-core  interaction,  Tliomas-Fermi,  and  multipole  expansion  atomic  force  con- 
stants. 
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